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reconstructor  to  handle  the  effects  of  intensity  dropout,  branch  points,  and  branch  cuts  is 
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I.  INTRODUCTION 


A.  MOTIVATION  AND  BACKGROUND 

Future  technologies  are  beginning  to  use  both  directed  energy  and  highly  accurate 
imaging  for  more  advanced  and  efficient  means  of  accomplishing  tasks,  which  previously 
were  thought  to  be  impossible.  The  use  of  high-power,  highly  accurate  optics  has  evolved 
with  both  space-based  and  ground-based  imaging  systems,  allowing  for  exploration  and 
discovery  into  areas  that  were  previously  untouched.  The  use  of  laser  energy  is  evolving 
to  allow  for  accurate  and  low  collateral  weapons,  in  addition  to  high-bandwidth 
communications  in  free-space  optical  communications.  The  initial  use  of  these  laser 
systems  was  subjected  to  the  environment  of  space,  which  had  its  own  non-atmospheric 
disturbances  and  issues.  Now,  that  technology  is  being  adapted  to  ground-based  systems 
which  experience  not  only  many  of  the  same  issues  space-based  encountered,  but  also  the 
new  difficulties  involved  with  energy  and  optical  propagation  through  an  enviromnent. 
The  possibilities  of  these  technologies  has  exploded  an  interest  amongst  scientists  and 
engineers  into  advancing  current  beam  control  and  high-precision  optical  capabilities. 

In  recent  years  the  Office  of  Naval  Research  (ONR)  has  undertaken  the 
challenges  posed  by  directed  energy  (DE)  weapons  with  the  creation  of  the  Directed 
Energy  Weapon  Program.  This  program  is  developing  a  directed  energy  laser  weapon 
that  will  increase  the  Navy’s  effectiveness  at  shooting  down  enemy  weapons  and/or 
hostile  craft  while  minimizing  collateral  damage.  The  Navy  has  identified  the  key 
challenges  of  a  DE  weapon  as  platform  stabilization,  imaging  and  correction  of  the  beam 
due  to  deep  atmospheric  turbulence,  and  simplifying  the  beam  control  system  operational 
cost  and  complexity.  Beam  control  for  directed  energy  weapons,  like  any  weapon,  is 
necessary  to  correctly  aim  and  fire  the  weapon  in  a  combat  environment.  A  primary 
difference  between  conventional  weapons  and  a  directed  energy  weapon  is  the  necessity 
to  maintain  a  highly  precise  and  power-dense  beam,  to  maximize  efficiency  and  lower 
dwell  time.  The  longer  the  beam  travels  through  the  atmosphere,  the  more  the  beam  is 
absorbed  and  scattered.  The  government  and  civilian  sectors  have  also  started  to  take  on 
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the  problems  involved  in  high-precision  imaging  and  optical-beam  control.  The  problems 
are  extremely  evident  in  the  programs  such  as  Hubble  and  James- Webb  space  telescopes. 

Precision  optical  beam  control  systems  can  be  broken  down  into  two  primary 
areas  of  control  and  correction;  jitter  control  and  adaptive  optics.  Jitter  control,  primarily 
a  concern  with  directed  energy  and  laser  communications,  focuses  on  the  beam  deviation 
induced  by  mechanical  and  atmospheric  vibrations  and  inaccuracies.  Adaptive  optics  was 
originally  used  by  ground  telescoped  for  correcting  image  distortion  due  to  the 
atmosphere.  It  is  also  used  for  correction  of  high  energy  laser  beam  aberrations  due  to  the 
atmosphere,  like  in  the  Airborne  Laser  (ABL).  Adaptive  optics  is  also  considered  to 
correct  beam  aberrations  due  to  imperfect  optical  surfaces. 

B.  THESIS  OBJECTIVE 

The  focus  of  this  research  is  to  investigate  new  methods  for  wavefront  estimation 
and  to  compare  with  other  commonly  used  methods.  Methods  and  tested  are  entirely 
implemented  in  a  simulation  based  environment.  The  research  focuses  on  the  ability  of 
these  different  wavefront  estimation  methods  to  deal  with  the  specific  issues  associated 
with  a  maritime  environment.  These  specific  issues  include  intensity  dropouts,  branch 
points,  and  branch  cuts.  The  images  from  each  reconstructor  are  compared  to  an  under 
sampled  version  of  the  original,  and  an  intensity-weighted  average  for  the  overall 
accuracy  of  the  reconstruction  is  calculated.  This  perfonnance  measurement,  commonly 
referred  to  as  a  Strehl  ratio,  is  used  as  the  primary  source  of  comparison  and  evaluation  of 
each  reconstructor.  Information  gained  from  this  research,  will  allow  for  confirmation  on 
whether  more  advanced  wavefront  estimation  methods  are  necessary  for  handling  the 
problems  associated  with  this  new  operating  environment. 

C.  THESIS  OVERVIEW 

Chapter  II  provides  background  information  and  the  basic  goals  of  an  adaptive 
optics  system,  as  well  as  the  basic  control  topology  associated  there  in. 

Chapter  III  provides  more  infonnation  on  the  issues  concerned  with  DE  and 
optical  imaging  propagation  through  the  atmosphere.  The  specific  problems  created  in  a 
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maritime  environment  are  introduced  the  principal  concepts  for  adaptive  optics  correction 
is  formulated. 

Chapter  IV  introduces  the  concepts  and  importance  associated  with  wavefront 
estimation,  and  its  role  in  the  adaptive  optics  control  process.  The  primary  means  of 
wavefront  sensing  is  introduced  and  modeled.  Then  the  three  wavefront  reconstruction 
and  estimation  methods  are  introduced  and  explained  including  the  two  commonly  used 
in  practices  today,  zonal  and  modal  estimation,  and  the  new  Noise  Variance  Weighted 
Complex  Exponential  (NVWCER)  reconstructor. 

Chapter  V  provides  simulated  results,  analysis,  and  comparison  of  all  three 
reconstruction  methods.  Their  ability  to  handle  the  specific  issues  that  lay  within  a 
maritime  environment  is  analyzed. 

Chapter  VI  provides  an  overall  summary,  conclusion,  and  recommendations  for 
possible  future  research. 
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II.  ATMOSPHERIC  TURBULENCE 


The  first  step  in  solving  a  problem  is  fully  understanding  the  problem  itself.  This 
section  presents  some  background  information  and  some  commonly  used  parameters 
used  to  describe  the  effects  of  turbulence.  There  are  two  primary  causes  for  turbulence  in 
the  atmosphere  as  a  laser  propagates  through.  The  primary  source  is  through  the  heating 
and  cooling  of  the  surface  of  the  earth,  which  in  turn  heats  and  cools  the  air  which 
changes  its  index  of  refraction.  These  differences  in  atmospheric  temperature  and  density, 
producing  differing  indexes  of  refraction,  cause  the  beam  to  diffract  disproportionately 
across  the  wavefront.  A  secondary,  and  less  common  source,  is  the  actual  heating  of  the 
air  by  the  laser  its  self,  in  turn  also  changing  the  index  of  refraction.  Changes  in  the  index 
of  refraction  cause  the  beam  to  diffract  and  dissipate  as  it  propagates  to  the  target, 
ultimately  creating  a  less  accurate  and  less  power-dense  beam.  AO  systems  focus  on 
correcting  the  disturbances  caused  by  the  first  source,  and  in  order  to  control  the  beam  it 
is  necessary  to  measure  the  beam  profile,  or  the  effects  of  the  turbulence  on  the  beam. 

A.  PARAMETERS  AND  MEASURMENTS 

There  are  two  main  paths  for  propagation  of  a  DE  beam:  vertical,  having  a  very 
high  slant  angle,  and  horizontal,  having  a  very  low  slant  angle.  Atmospheric  turbulence 
for  vertical  paths,  i.e.,  ground-to-space  architectures,  has  been  modeled  and  simulated 
using  Kolmogorov  turbulence  theory  and  statistics. 

Kolmogorov  theory  assumes  that  all  small  scale  turbulent  motions  are  directly 
related  to  parent  large  scale  turbulent  motions,  and  are  statistically  homogenous  and 
isotropic  (Roggemann  &  Welsh,  1996)  This  idea  comes  from  the  assumption  that  a 
turbulent  atmosphere  is  made  of  various  eddy  sizes,  but  the  larger  eddies  continuously 
break  down  into  smaller  eddies,  unifonnly  distributing  the  energy  until  the  viscosity  of 
the  fluid  allows  for  complete  dissipation  of  all  kinetic  energy  to  internal  energy. 
Kolmogorov  is  then  able  to  mathematically  describe  the  spatial  frequency  characteristic 
of  the  index  of  refraction  over  a  propagation  path  (Corley,  2010).  The  scale  of  those  eddy 
sizes  which  Kolmogorov  covers  is  referred  to  as  the  inertial  subrange.  Any  eddy  larger 


5 


then  L0  is  assumed  not  to  be  homogenous  in  the  atmosphere  and  any  eddy  smaller  then 

Z0  dissipates  the  kinetic  energy  as  heat  instead  of  transferring  it  to  smaller  eddy.  Using  the 

11 

inertial  subrange  where  —  «  k  «  —  the  power  spectral  density  (PSD)  of  the  index  of 

^0  ^0 


refraction  in  air  is 


<Pn(k,z)  =  0.033  C^(z)/c" 


(2n\  2n 

—  I,  is  related  to  the  isotropic  scale,  /,  by  l  —  — ,  z  is  the 

distance  from  the  aperture,  and  C%  is  a  measure  of  the  strength  of  the  turbulence 

_2 

(Andrews,  2004).  Weak  turbulence  is  associated  with  a  C%  of  10_17m_3  or  less,  while 

_2 

strong  turbulence  generally  has  a  C%  value  of  10_13m_3  or  greater.  One  issue  though  is 
that  values  vary  with  height  above  ground,  so  it  is  necessary  to  have  altitude 
dependent  model  or  profile.  These  profiles  do  exist,  and  have  been  created  from 
extensive  data  collection  and  is  based  on  geographic  location,  temperature,  particle  count, 
and  many  other  factors.  One  concern  however  is  that  these  profiles  have  primarily 
focused  on  vertical  paths  due  to  previous  mission  requirements,  so  new  data  is  being 
collected,  with  a  focus  on  the  maritime  environment  and  lower  slant  angles  (Corley, 
2010). 


Another  common  parameter  used  to  describe  atmospheric  turbulence  is  Fried’s 
parameter,  or  atmospheric  coherence  length,  r0.  Fried’s  parameter  describes  a  seeing 
condition  at  a  given  distance,  stating  that  a  telescopes  resolution  is  diffraction  limited  at 
any  diameter  larger  thanr0  (Andrews,  2004).  Fried’s  parameter  is  measured  in 
centimeters,  and  is  formulated  as 

0.42 sec(Qk2  f  C^{z)dz 
Jo 

where  C,  is  the  zenith  angle,  L  is  the  distance  from  the  source  to  telescope  aperture  along 
the  primary  axis,  and  k  is  the  wavenumber.  Values  of  5  cm  or  less  represent  strong 
turbulence,  where  any  values  over  25  cm  represent  very  weak  turbulence,  and  good 
seeing  conditions  (Wilcox,  2009) 
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Kolmogorov,  Fried,  and  C%  adequately  describe  the  spatial  structure  of  the 
atmosphere  and  its  turbulence,  however  it  does  not  account  for  the  time  varying  nature  of 
the  atmosphere  experienced  with  long  dwell  times  or  slewing  engagements.  While  there 
are  many  ways  to  simulate  or  characterize  the  time  varying  nature  of  atmospheric 
turbulence,  the  Greenwood  time  constant  is  a  good  measure  of  the  temporal  coherency  of 
the  atmosphere.  Greenwood  measures  a  time  constant  over  which  a  propagation  path  can 
be  assumes  to  have  a  constant  index  of  refraction,  before  needing  to  be  recalculated.  The 
time  constant  is  measured  by 

3  i 

cos 5 (£)  0.314r0  (D\c 

T° =  r  ~  5  1 

2.91/c2  C%(K)VHh)dh 

where  V(h)  is  the  wind  velocity  as  a  function  of  altitude,  but  can  also  represent  the 
slewing  rate  of  an  moving  engagement  scenario. 

In  order  to  analyze  the  effects  atmospheric  turbulence  is  having  on  a  beam  or 
wavefront  it  is  useful  to  generate  a  mathematical  representation,  rather  than  just  a  discrete 
set  of  values.  Since  it  is  assumed  a  wavefront  is  generally  continuous  over  a  two 
dimensional  surface  wavefront  for  certain  atmospheric  conditions,  it  makes  sense  to 
express  the  wavefront  as  a  linear  combination  of  polynomial  based  functions  which  are 
orthonormal  over  a  an  entire  optical  surface.  The  set  of  linear  polynomials  used  in  this 
research  are  Zernike  polynomial  represented  by 

m 

Phase(p,9 )  =  I  aiZi  (p,  0) 

i= 1 

where  m  is  the  number  of  desired  modes  and  ai  is  the  weight  associated  with  each  mode 
(Allen,  2007).  The  accuracy  of  a  Zernike  polynomial  representation  is  a  function  of  both 
the  amount  of  disturbance  in  the  wavefront  and  how  many  modes  are  chose  for  the 
Zernike  polynomial  to  represent  the  disturbance.  Theoretically  as  m  approaches  infinity 
you  can  model  any  continuous  wavefront,  however  is  beneficial  to  choose  the  correct 
number  of  modes  based  on  two  criteria:  how  turbulent  is  the  data  itself  and  how  much 
accuracy  can  your  sensors  or  controllers  actually  use. 


7 


B. 


MARITIME  ENVIRONMENT  AND  DEEP  TURBULENCE 


Moving  to  the  maritime  environment  creates  compounded  errors  when  compared 
the  traditional  ground-to-space  operating  procedures.  First  is  the  effect  of  a  significantly 
longer  propagation  distance  through  the  atmosphere.  Secondly,  models  already  created 
based  on  C%  measurements  were  for  vertical,  versus  horizontal,  paths.  These  problems 
combined  lead  to  a  need  for  more  robust  adaptive  optics  control. 

As  previously  discussed  by  Melissa  Corely  (Corely,  2010),  there  is  work  currently 
underway  to  create  the  new  atmospheric  disturbance  profiles,  but  there  is  no  model  for 
horizontal  paths  that  matches  the  familiarity  and  accuracy  of  the  Kolmogorov  model  for 
vertical  paths  (Corely,  2010).  Agencies  and  organizations  like  SPAWAR  in  San  Diego, 
NPS  in  Monterey,  the  Naval  Research  Lab  (NRL),  University  of  Puerto  Rico,  University 
of  Florida,  and  Michigan  Tech  have  all  begun  work  to  develop  and  characterize 
parameters  like  C^and  ro  for  horizontal  paths  over  bodies  of  water.  NPS  has  even 
developed  an  initial  model,  the  Navy  Surface  Layer  Optical  Turbulence  (NSLOT)  model 
which  models  as  a  function  of  local  air  and  sea  temperatures. 

While  work  is  being  done  to  create  accurate  models,  that  is  only  one  of  the  many 
issues  with  maritime  environment  beam  propagation.  The  presence  of  higher  amounts  of 
aerosols,  large  humidity  and  temperature  fluctuations,  and  wave  motion  all  contribute  to 
higher  levels  of  scintillation.  Scintillation  is  intensity  fluctuations  or  random  changes  in 
the  index  of  refraction,  and  this  causes  a  large  problem  with  adaptive  optics.  One  it 
degrades  the  quality  and  energy  density  on  target,  but  it  also  makes  it  difficult  for  sensors 
to  accurately  estimate  the  wavefront  and  use  it  for  correction. 

C.  BRANCH  POINTS  AND  BRANCH  CUTS 

The  scintillation  effects  and  intensity  dropouts,  experienced  in  the  presence  of 
deep  turbulence,  create  nulls  in  the  intensity  making  it  impossible  for  a  wavefront  sensor 
to  obtain  a  reading  at  those  points.  These  null  readings  in  intensity  are  also  shown  to  lead 
to  a  non-singular  valued  function  when  calculating  the  phase.  These  non-singular  value 
functions  mean  when  a  closed-loop  path  around  a  subjective  point  is  calculated,  one  gets 
a  non-zero  value,  which  is  also  known  as  a  branch  point.  Branch  points  have  either  a  net 
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negative  or  positive  discontinuity,  and  negative  branch  points  must  be  connected  to  a 
corresponding  positive  branch  point,  through  a  2n  discontinuity  known  as  a  branch  cut. 

Branch  points  and  branch  cuts  drastically  decrease  the  efficiency  with 
which  a  traditional  AO  system  can  measure  and  compensate.  That  is  where  the  focus  of 
this  research  will  be;  attempting  to  accurately  handle  and  model  branch  points  and 
intensity  dropouts  in  wavefront  estimation  for  an  AO  system.  Classical  wavefront 
reconstructor  use  a  least  squares  fit  method,  which  in  a  sense  ignores  the  hidden  phase 
and  glances  over  branch  points  and  branch  cuts.  A  least  mean  squares  reconstruction 
method  can  be  viewed  as  forming  an  estimate  of  a  phase  point  through  adding  up  the 
phase  differences  along  every  possible  path  from  a  reference  point  to  the  desired  point 
and  averaging  them  together  (Fried,  1998).  This  method  works  well  in  a  noise- weighted 
sense,  however  in  the  presence  of  branch  points  different  paths  can  lead  to  different 
values,  both  of  which  can  be  equally  correct.  The  issue  is  that  branch  points  in  least 
squares  estimates  of  the  wavefront  have  a  degrading  effect  over  a  wide  area  of  the 
aperture  surrounding  the  branch  points,  as  multiple  paths  will  go  through  this  region. 
Figure  1  shows  the  effects  of  branch  points  and  branch  cuts  on  a  path  following 
understanding  of  least  mean  squares  estimators. 
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Figure  1.  Phase  calculation  paths  leading  to  different  phase  values  (From  Pellizzari 

2010) 
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While  it  may  seem  incorrect  at  first  to  arrive  at  two  different  phase  values  for  the 
same  point,  they  are  both  equally  correct.  The  error  stems  from  an  assumption  made  by 
least  squares  estimators,  and  that  is  that  the  measured  phase  differences  are  purely  a 
function  of  a  gradient  of  phase  function.  The  measured  phase  differences  need  to  be 
treated  similarly  as  an  electromagnetic  field  with  flowing  electricity,  and  that  is  as  a  sum 
of  a  scalar  potential  plus  the  curl  of  a  vector  potential  (Fried,  1998).  This  view  of  the 
phase  gradient  as  containing  both  a  real  and  imaginary  part  goes  back  to  the  governing 
equations  for  electromagnetic  propagation  known  as  Maxwell’s  equations.  A  new 
reconstruction  method  is  explored,  which  takes  into  account  the  hidden  phase  aspect,  and 
identifies  the  existence  and  location  of  branch  points  and  branch  cuts  in  its  wavefront 
calculations.  When  branch  points  occur,  it  is  ideal  for  them  to  be  located  in  areas  of  low 
intensity  of  the  aperture  to  minimize  the  effects. 


10 


III.  ADAPTIVE  OPTICS 


Adaptive  optics  technology  is  a  multidisciplinary  field  which  combines  expertise 
from  engineers  and  scientist  in  wide  range  of  fields  including  mechanics,  optics, 
materials  science,  chemistry,  physics,  and  control  theory  (Allen,  2007).  The  primary 
theories  associated  with  adaptive  optics  have  not  changed  over  the  years;  however  they 
are  becoming  more  refined  and  focused.  Initially  the  components  used  in  adaptive  optics 
like  cameras,  defonnable  mirrors,  and  computing  capabilities  were  very  coarse  by 
today’s  standards  and  only  a  certain  level  of  control  precision  could  be  detected  and 
applied.  With  large  advances  in  this  same  technologies  allowing  for  greater  precision  and 
detail,  compounded  with  the  introduction  of  new  technologies  like  segmented  mirrors, 
the  requirements  for  advanced  adaptive  optics  systems  is  growing  rigorously.  New 
developments  in  technology  and  theory  are  leading  scientists  and  engineers  to  find  new 
applications  for  AO,  which  in  turn  creates  even  more  specific  demands  for  an  AO 
system;  revolving  in  a  continuous  loop. 

The  goal  of  an  AO  system  is  to  actively  improve  the  capability  of  an  optical 
system  in  the  presence  of  imperfect  operating  conditions.  An  AO  system  continuously 
attempts  to  actively  compensate  for  two  main  sources  of  optical  aberrations: 
imperfections  in  the  surface  of  a  lens  in  the  optical  system  and  wave  front  distortion 
produced  by  Earth’s  atmosphere.  AO  was  originally  used  for  space  imaging  systems  to 
handle  the  aberrations  induced  by  Earth’s  atmosphere.  Now  AO  is  being  employed  in  DE 
energy  applications  for  both  weapons  and  communications,  in  which  coherency  of  the 
wavefront  phase  is  necessary  for  precision  and  energy  density. 

A.  TYPICAL  ADAPTIVE  OPTICS  SYSTEMS 

Most  adaptive  optics  system  follows  a  setup  based  on  the  schematic  show  in 
Figure  2,  which  consist  of  three  components.  The  first  is  a  wavefront  sensor,  which 
provides  some  measurement  of  the  optical  wavefront;  that  can  either  be  discrete  slope 
values  or  phase  values  at  given  points.  The  second  component  is  the  computer  controller 
which  attempts  to  reconstruct,  or  estimate,  the  phase  of  the  wavefront  and  compute  a 
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control  signal  for  the  deformable  mirror.  The  third  component  is  a  deformable  mirror, 
which  a  membrane  like  optical  surface  mounted  on  tiny  electrical-actuators  which  can  be 
used  to  defonn  the  mirror  to  change  the  phase  of  the  reflecting  wavefront.  The  wavefront 
sensor  and  computer  controller  work  together  to,  as  accurately  as  possible,  estimate  what 
the  wavefront  is  and  figure  out  which  points  need  to  be  either  delayed  or  advanced  so  that 
the  wavefront  is  as  flat  as  possible.  This  research  will  be  focusing  on  the  sensor  and 
computer  controller  portion  of  the  AO  system  and  the  specifics  of  the  sensor  will  be 
introduced  in  the  following  chapter. 

It  is  important  to  note  that  the  reference  beam  can  come  from  a  wide  array  of 
sources;  however  for  astronomical  applications  it  is  common  to  use  a  known  star.  A  star 
imaged  from  an  extremely  far  distance  allows  it  to  appear  as  a  point  source,  so  it  is  easy 
to  determine  if  any  blurring  or  aberrations  are  present,  and  allow  for  correction. 


tt 


Figure  2.  Typical  AO  System  Setup  (From  Allen  2007) 


B.  BASIC  ADAPTIVE  OPTICS  CONTROL  CONCEPTS 

While  this  research  focuses  on  the  wavefront  reconstruct  portion  of  AO  system,  in 
order  to  have  a  full  understanding  of  the  motivation  behind  the  problem,  it  requires  a 
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basic  understand  of  the  control  principals  associated  with  it.  While  the  control  computer 
works  with  the  wavefront  sensor  to  determine  the  errors,  in  order  for  the  computer  to 
compute  a  control  signal  it  needs  some  sort  of  relationship  between  the  deformable 
mirror  and  the  wavefront.  If  one  considers  a  classical  negative  feedback  system 
illustrated  in  Figure  3  then  the  wavefront  sensor  can  be  interpreted  as  a  state  estimator, 
the  optical  wavefront  is  the  “plant”  that  is  controlled,  and  the  deformable  mirror  is  the 
actual  control  element  (Tyson,  2004).  In  Figure  2  there  is  also  a  reconstructor  block,  and 
this  is  what  provides  that  relationship  between  the  deformable  mirror  and  the  wavefront 
sensor. 


Tuibulence 


m  *-  > 


Reconstructor 
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Figure  3.  Classical  AO  Control  System 


In  a  simple  overview  a  wavefront  estimator,  which  can  a  multitude  of  different 
sensors,  represents  the  wavefront  as  a  vector  of  discrete  values.  These  values  can  either 
be  local  slope  measurements,  local  phases,  or  even  coefficients  associated  with  equations 
like  Zernike  polynomials.  Recalling  that  a  deformable  mirror  has  a  discrete  number  of 
actuators,  either  electrostrictive  (PZT)  or  magnetostrictive  (PMN),  a  few  other 
assumptions  need  to  be  made  as  well.  In  modern  AO  systems  it  is  assumed  there  is  static 
coupling  between  actuators,  however  there  is  not  dynamic  coupling  in  the  deformable 
mirror  (DM).  The  response  time  of  a  DM  is  so  quick  the  relationship  between  voltage  and 
actuator  movement  is  treated  as  a  step  response.  It  is  also  important  to  note  that  the 
relationship  between  voltage  and  actuator  movement  can  be  considered  linear  with  a 
small  displacement  range.  With  these  assumptions  the  concept  of  a  poke,  or  influence, 
matrix  can  be  constructed. 

The  deformable  mirror  has  a  discrete  amount  of  actuators  and  the  wavefront 
sensor  outputs  a  discrete  amount  of  coefficients  or  slopes/phases,  so  each  actuator  is  fully 
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triggered  and  the  corresponding  change  to  the  outputs  of  the  wavefront  sensor  are  logged. 
This  is  done  sequentially  and  individually  for  each  actuator,  while  holding  all  others  at 
zero.  These  measurements  are  then  used  to  create  a  poke  matrix  whose  columns 
correspond  to  wavefront  sensor  changes  associated  with  that  actuator,  and  the  rows 
correspond  to  the  same  coefficient  or  slop  reading  across  each  actuator.  The  structure  of  a 
typical  poke  matrix  using  a  Shack-Hartmann  wavefront  sensor  (SHWFS),  which  outputs 
local  x  and  y  sloped  measurements  at  discrete  points,  is: 


-Xu 

...  xlk- 

X21 

...  X2k 

xnl 

xnk 

Til 

-  y-ik 

721 

—  72  k 

-7n  1 

7nk- 

where  n  is  the  number  of  discrete  x  and  y  slope  measurements  and  k  is  the  number  of 
actuators  on  the  deformable  mirror.  In  the  case  of  a  sensor  with  n  slope  measurement 
points,  and  k  actuators  the  matrix  is  2n  x  k.  Using  the  influence  matrix  the  relationship 
between  command  voltage,  c,  and  slope  output,  s,  from  the  senor  is  defined  as 

s  =  I'c 

Given  the  control  architecture  of  AO  system,  usually  the  slope  error  is  known, 
and  a  actuator  command  voltage  is  desired.  Normally  the  multiple-input-multiple-output 
(MIMO)  AO  system  has  more  inputs  (sensor  readings)  than  outputs  (actuator 
commands),  and  is  therefore  over  detennined  and  creating  a  non-square  matrix.  Taking 
the  pseudo  inverse  of  the  poke  matrix,  F+,  and  pre-multiplying  the  slope  error 
measurements  creates  a  least  squares  solution  to  the  minimization  of  the  wavefront  error 
correction: 

r+s  =  r+rc  =  c 

In  Figure  3  the  controller  implemented  is  a  simple  integrator,  which  uses  a  gain, 
g,  to  iteratively  calculate  a  new  control  using  the  following  law 

Cnew  —  Cold.  9^  S 

This  is  just  one  example  of  a  simple  controller,  but  any  control  method  can  be 
implemented. 
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IV.  WAVEFRONT  ESTIMATION 


The  ability  to  reconstruct  and  estimate  the  wavefront  is  essential  for  a  successful 
AO  system.  In  traditional  wavefront  sensing  for  optical  purposes,  the  phase  determination 
can  take  anywhere  from  hours  to  days,  however  adaptive  AO  systems  operate  at  a  higher 
demand.  AO  systems  are  used  for  real-time  corrections  and  need  to  operate  at  rates  up  to 
several  hundred  of  hertz  and  spatial  resolution  as  fine  as  1/100  of  the  aperture  diameter. 
There  are  two  methods  for  sensing  a  wavefront;  direct  and  indirect  approach  (Tyson, 
1998).  The  direct  approach  takes  sensor  data  and  attempts  to  explicitly  solve  for  the 
phase  of  the  wavefront  which  is  used  as  information  for  correcting  unwanted 
components.  The  indirect  method  directly  relates  the  sensor  data  to  a  control  signal, 
bypassing  the  need  to  explicitly  solve  for  the  phase  (Tyson,  1998).  While  the  indirect 
may  be  quicker,  the  direct  method  provides  greater  detail  of  the  errors  in  the  wavefront, 
therefor  determination  of  the  method  used  is  completely  situation  dependent. 

A.  WAVEFRONT  SENSOR 

There  are  various  sensors  used  for  wavefront  detection;  this  research  makes  use  of 
a  Shack-Hartmann  wavefront  sensor  (SHWFS)  design.  A  SHWFS  consist  of  two 
dimensional  array  of  lenslet  placed  in  a  grid  pattern  in  front  of  an  imaging  sensor,  as 
shown  in  Figure  3. 


Each  lenslet  acts  as  an  aperture,  causing  the  incoming  light  to  converge  into  spots 
on  the  CMOS  sensor.  The  location  of  each  spot  is  directly  proportional  to  the  average 
slope  of  the  wavefront  passing  through  the  corresponding  lenslet  (Primot,  2003).  A 
perfectly  coherent  wavefront  is  used  to  produce  a  reference  image  for  centroid  locations. 
The  deviation  of  each  centroid  for  an  aberrated  wavefront  from  each  reference  location  in 


the  x  and  y  directions  is  directly  proportional  to  the  local  x  and  y  slopes,  a  and  P, 
respectively.  Figure  4  shows  the  relationship  between  the  slope  of  the  wavefront  and  the 
displacement  of  the  centroid.  The  displacement  of  the  centroid  is  related  to  the  slope  by: 

2n 

aiJ  = 

2n 

hi  = 

where  lambda  is  the  wavelength  of  light  passing  through,  f  is  the  focal  length  of  the 
lenslet  array,  and  delta  x  and  delta  y  are  the  shifts  of  the  centroids  in  each  dimension, 
respectively  (Zhu,  Sun,  Bartsch,  Freeman  &  Fainman,  1999). 


The  centroids  are  calculated  for  each  array  of  pixels  corresponding  to  the 
equivalent  sub-lenslet.  The  Centroid  calculation  is  done  through  a  summation  of  the  pixel 
location  from  the  given  axis  weighted  by  the  corresponding  pixel  intensity,  defined  as 


Zu=1Z”=i/(m,  v) 
ZSUZSLi  v*i(u,v) 
ESLiZSUK  u,v) 


where  i  and  j  specify  the  given  lenslet  and  u  and  v  correspond  to  the  local  pixel  indexes 


of  the  given  supabperature  with  (1,1)  being  the  bottom  left  corner. 


The  algorithm  is  carried  out  on  a  small  section  of  pixels  from  the  CMOS  sensor, 
corresponding  to  the  given  lenslet.  It  is  assumed  the  very  small  gaps  between  the  circular 
lenslet  arrays  create  minimal  noise  on  the  CMOS  sensor  through  leakage,  therefore  it  can 
be  ignored.  With  the  corresponding  slope  measurements  then  either  a  control  signal  can 
be  calculated  for  the  DM,  using  the  indirect  control  method,  or  the  phases  of  the 
wavefront  can  be  reconstructed  if  a  direct  approach  is  used. 
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B. 


SHWFS  OUTPUT  SIMULATION  METHOD 


The  research  was  conducted  in  an  entirely  simulated  environment,  so  a  method 
for  simulating  the  SHWFS  is  needed.  Using  WAVEPROP,  developed  by  the  Optical 
Sciences  Company  (tOSC),  to  simulate  the  effects  of  a  turbulent  environment  on  a  beam, 
through  the  use  of  Kolmogorov  statistics,  an  incident  phase  and  amplitude  screen  is 
created.  These  phase  and  amplitude  screens  simulate  what  will  be  seen  at  the  face  of  a 
SHWFS.  The  SHWFS  is  a  planar  array  of  multiple  lenses  with  a  sensor  placed  at  the 
focal  point;  electro-optical  theory  can  be  used  to  simulate  the  effects  of  each  lenslet.  In 
wave  optics  each  lens  can  be  modeled  with  a  mathematical  expression  of  a  phase  device 
(Schmidt,  2010) 


lensix.y )  =  w(x,y )  exp 


where  k  is  the  wave  number,  /  is  the  focal  length  and  w(x,y)  represents  the  aperture 
function 


w(x,y)  = /(x)  =  j1'  1*1  <  2  'lyl<2 

lO,  other 

where  d  is  the  size  of  the  lens.  Fourier  optics  provides  a  simple  calculation  to  evaluate  the 
intensity  distribution  at  the  focal  plane  as 


l{xf,yf)a  j!  i/t(x,  y)w(x,y  )  exp  -i  (xxf  +  yyf)  dxdy 
=  |F7’{T/;1(x,y)w(x,y)}|2 

where  vp  is  the  complex  wave  distribution  of  the  phase  and  amplitude  screen  at  the  face  of 
the  lenslet  array  defined  as 


xp(x,y)  —  A(x,y~)  exp(i  *  0(x,y)) 

With  A(x,y)  being  the  amplitude  and  9(x,y)  being  the  phase  screens,  respectively.  While 
Fourier  optics  provide  a  method  of  simulating  the  optical  effects  of  each  lens,  it  has  been 
proven  that  the  equivalent  slope  output  through  a  centroiding  algorithm  of  a  SHWFS  is 
directly  proportional  to  the  intensity  weighted  average  of  the  gradient  between  each  pixel 
across  the  entire  sensor  (Fried  1997).  This  relates  to  the  centroiding  algorithm  used  to 
calculate  the  displacement  of  a  SHWFS,  in  which  areas  of  higher  intensity  have  a  greater 
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effect  on  the  centroid  location.  The  equivalent  x  and  y  gradient,  or  tilt,  can  be  calculated 
for  each  supabperature  as 

771-1 


v — 1  v — 1  m 

sx(uj )  =  2^ 


U-l 

771 


\  1  ^ — 1 771—  1 

Syd.i)  =  X  ^ 


( arg(xp(u  +  1,v)i/j*(u,  i;)))  ( abs(xp(u  +  1,  v)ip*(u,  v))j 

Z«=i  ZSLi  ( abs(xp(u  +  1,  v)xp*(u,  v))) 
[arg{xp{u,v  +  1  ')rp*(u,v')'fj  ( abs(i/j(u,v  +  l)i/F(it,  v))) 


u=i  '  *  Xu=i'Z™=i(abs(xP(u’v  +  i^CzmO)) 

with  u  and  v  corresponding  to  the  indices  of  the  phase  and  amplitude  screens  which 
correspond  to  the  area  covered  by  the  given  lenslet  and 


( Im(x)\ 
arg(x)  =  atan(j^) 


A  simple  test  algorithm  was  run  to  ensure  the  proper  functioning  of  the  SHWFS 
simulator.  A  phase  plane  shown  in  Figure  5  is  defined  by  the  equation 

0(x,y)  —  x2  +  y3 

which  generated  the  following  phase  screen 


o  o 


Figure  5.  SFIWFS  Slope  Simulate  Test  Phase  Screen  (z=x2+y3) 

The  weighted  averaging  was  used  to  simulate  the  slopes  generated  by  a  64x64 
lenslet  SFIWFS,  and  Figure  6  shows  the  results. 
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0.3 


X  Slope  Data  Y  Slope  Data 


Figure  6.  SFIWFS  Output  Simulation  for  X  and  Y  slopes  based  on  Fig.  5  Phase 


From  basic  mathematical  differentiation  of  the  phase  equation,  it  makes  sense  that  the  x 
slope  follows  a  linear  path  and  y  data  follows  a  quadratic  (second  order)  growth  path.  A 
normally  distributed,  random  amplitude  screen  was  generated  and  applied  for  the 
algorithm. 

C.  ESTIMATION  (RECONSTRUCTION)  METHODS 

As  discussed  in  section  III  there  are  two  methods  of  actually  implementing  the 
slope  data  from  a  SFIWFS  for  control  of  a  DM,  either  direct  or  indirect,  and  wavefront 
estimation  is  used  with  the  direct  method.  There  are  numerous  methods  to  actually  go 
from  slope  data  to  an  estimated  phase.  Two  common  methods,  typically  used  in  low 
turbulence  environment,  are  zonal  and  modal  estimation  (Baker,  2005).  Zonal  attempts  to 
estimate  discrete  phase  points  in  each  zone,  while  modal  attempts  to  solve  for  the 
coefficients  of  polynomial  expansion  like  functions.  A  third  method,  the  focus  of  this 
research,  is  a  noise-variance  weighted  complex  exponential  reconstructor,  which  is  a 
recursive  noise-weighted  averaging  algorithm  which  use  complex  phase  differentials 
instead  of  just  slopes. 

First  the  different  geometries  need  to  be  introduced  which  relate  the  slope 
measurements  of  the  SHFWS  to  the  calculated  phase  points.  Figure  7  shows  the  differing 
geometries  where  the  lines  indicated  the  slope  or  phase  difference  values  and  the  dots 
represent  the  calculated  phase  values.  Zonal  and  Modal  reconstruction  methods  operate 
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on  Southwell  geometry,  while  the  NVWCER  operates  on  a  combined  Hudgin  and  Fried 
geometry. 


(A) 
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(B) 

1111 
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4  4  4 

4  4  4  4 

•  -  •  -  •  -  • 

•  •  •  • 

Figure  7.  Common  Reconstruction  Geometries  (a)  Southwell  (b)  Hudgin  (c)  Fried 

(From  Southwell  1980) 


1.  Zonal  Reconstruction 


Using  a  SHWFS  produces  incident  x  and  y  slopes  at  a  single  aperture,  so  the 

geometric  configuration  in  Figure  7(a)  is  used.  The  first  step  is  to  formulate  a  relationship 

between  the  phases  and  slopes.  If  quadratic  interpolation  is  used  to  define  the  dependence 

of  the  phase  in  the  x  and  y  direction  then  the  following  two  polynomials  are  created. 

0  =  c0  T  cxx  +  c2x2 
<f>  =  c0  +  c1y  +  c2y2 

with  relationship  for  the  dependency  of  the  phase  on  the  x  and  y  locations,  the 

differentiation  of  each  equation  should  provide  the  equivalent  dependency  equations  for 

the  x  and  y  slopes.  Taking  the  derivative  of  the  equations  with  respect  to  x  and  y  yields 

Sx  —  ci  +  2  c2x 
Sy  =  c1  +  2  c2y 

Since  two  slope  measurements  are  given  for  each  phase  point  it  allows  for  the 
determination  of  ci  and  C2.  A  relationship  between  adjacent  phases  and  slope 
measurements  is  created;  the  average  slope  of  two  adjacent  phase  points  is  equal  to  the 
difference  of  the  two  phase  points  divided  by  the  length  h  (Southwell,  1980). 


c*  .  +  c* 


2 


0t+lj  0ij 

h 


i  =  1:  (IV  -  1)  ;  j  —  1:N 
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5 


y 

i,j+ 1 


2 


4>i,j+ 1  <Pi,j 

h 


i  =  l:N;  j  =  1:(N  -  1) 


Each  lenslet  represents  an  equal  sub-area  of  the  entire  aperture,  h2,  so  /z=D/N 
where  D  is  the  diameter  of  the  entire  aperture  and  N  in  the  number  of  lenslets  along  one 
side.  In  order  to  simultaneously  solve  all  for  all  phase  points  in  a  least  squares  solution,  a 
matrix  interpretation  of  the  operation  is  created.  First  a  vector  of  the  slopes,  of  length 
2N“,  is  created  with  the  x  slopes  first,  then  y  slopes. 


2 

Next  a  sparsely  populated  matrix  D,  of  size  2N(N-1)  by  2N“,  is  created  which  performs 
the  adjacent  slope  averaging.  Then  a  matrix  A,  of  size  2N(N-1)  by  NT,  is  created  which 
calculates  the  adjacent  phase  differences.  Using  matrix  formulation  the  phase  and  slopes 
are  related  by 

[D]S=  [A\4> 

Given  that  D,  S,  and  A  are  known  the  phases  can  be  calculated  by  multiplying 
both  sides  by  the  inverse  of  A.  However,  since  there  are  more  slope  measurements  then 
phase  points  the  equation  is  over  determined  and  A  is  not  a  square  matrix,  so  the  pseudo 
inverse  must  be  taken,  which  results  in  a  least  squares  solution 

<? = ww 

where 

i/i]t «  nra.4H.4f)-1 


2,  Modal  Reconstruction 


In  modal  reconstruction  the  wavefront  is  described  with  coefficients  of  multiple 
different  optical  modes,  expressed  as  an  optical  wavefront  expansion;  similar  to  that  of  a 
polynomial  expansion.  It  makes  use  of  Southwell  geometry,  like  modal  reconstruction, 
assuming  coincident  slope  measurements.  There  are  numerous  mathematical  models  to 
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represent  the  different  modes,  and  this  thesis  focuses  on  the  use  of  Zernike  polynomials. 
Zernike  polynomials  have  two  key  characteristics:  they  are  orthogonal  over  a  unit  disk 
and  also  expandable  to  a  limitless  degree.  The  orthonormality  of  the  polynomials  means 
that  there  derivative  also  exists  over  the  entire  unit  disc,  which  becomes  useful  when 


attempting  to  fit  slope  data.  Zernike  polynomials  represent  a  wavefront  a  desired  m  and  n 
degree  through  (Tyson  and  Frazier,  2004) 


0(r,  9) 


^An0 K  Q  +  X  X  ( Anm  cos(md )  +  Bnm  sin(m6 ))  (£) 


n=  1  m= 1 


where 


n-m 

,r,  =  y  (-l)s((n-s)Q  ,r.  n-2s 

It  is  important  to  note  that  Zernike  polynomials  exist  in  polar  coordinates  and 


exist  over  a  unit  disk.  The  radius  is  normalized  with  the 


term,  where  r  is  the  radius  at 


a  given  point  and  R  is  the  radius  of  the  entire  aperture.  The  first  21  modes  are  listed  in 
Table  1. 
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# 

n 

m 

Polynomial 

0 

0 

0 

1 

1 

1 

1 

p  Cos [0] 

2 

1 

1 

p  Sin [0] 

3 

1 

0 

-  1  +  2  p2 

4 

2 

2 

p2  Cos  [2  0] 

5 

2 

2 

p2  Sin  [2  0] 

6 

2 

1 

p  (-2  +  3  p2)  Cos [0] 

7 

2 

1 

p  (-2  +  3  p2)  Sin [0] 

8 

2 

0 

1  -  6  p2  +  6  p4 

9 

3 

3 

p3  Cos [ 3  0] 

10 

3 

3 

p3  Sin [ 3  0] 

11 

3 

2 

p2  (-3  +  4  p2)  Cos  [2  0] 

12 

3 

2 

p2  (-3  +  4  p2)  Sin  [2  0] 

13 

3 

1 

p  (3  -  12  p2  +  10  p4)  Cos  [0] 

14 

3 

1 

p  (3  -  12  p2  +  10  p4)  Sin  [0] 

15 

3 

0 

-1  +  12  p2  -  30  p4  +  20  p6 

16 

4 

4 

p4  Cos [ 4  0] 

17 

4 

4 

p4  Sin [ 4  0] 

18 

4 

3 

p3  (-4  +  5  p-)  Cos  [3  0] 

19 

4 

3 

p3  (-4  +  5  p-)  Sin[30] 

20 

4 

2 

p2  (6  -  20  p2  +  15  p4)  Cos  [2 

Table  1 .  First  2 1  Zernike  Polynomial  Modes  (From  Wyant  2003) 

The  SHWFS  modeled  in  this  research  relies  on  a  square  grid  arrangement  which 
introduces  two  issues;  the  need  for  a  proper  coordinate  system  and  fitting  a  circle  to  the 
square.  Zemike  polynomials  can  be  transformed  into  Cartesian  coordinates  through  the 
use  of  two  relationships, 


-  yjx2  + 


r 


6  —  arctan  ) 


and  their  equivalent  equations  are  show  in  Table  2.  Similarly  to  polar  coordinates, 
Cartesian  coordinates  are  normalized  as  well  where  x  and  y  are  both  normalized  to  satisfy 

[p  =  yjx2  +  y2]  <  1 
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# 

n 

m 

Polynomial 

Term 

0 

0 

0 

1 

Piston 

1 

1 

1 

X 

X-Tilt 

2 

1 

1 

y 

Y-Tilt 

3 

1 

0 

-1  +  2^+/) 

Focus 

4 

2 

2 

x2-y2 

Astigmatism  plus 
defocus 

5 

2 

2 

2  xy 

Astigmatism  plus 
defocus 

6 

2 

1 

-2x  +  3x(x2  +y2) 

Coma 

7 

2 

1 

-2y  +  3y(x2+y2) 

Tilt 

8 

2 

0 

1_6(.Y2  +  j;2)  +  6(A.2+;);2j2 

Third-Order 

Spherical  and  Focus 

9 

3 

3 

x3  -  3  AY  2 

Fifth-Order 

Aberration 

10 

3 

3 

3  x2y-y3 

Fifth-Order 

Aberration 

11 

3 

2 

-3x2  +  3y2  +  4x2  (x2  +  y2 )  -  4y2  (x2  +  y2 ) 

Fifth-Order 

Aberration 

12 

3 

2 

-6ay  +  8xv  ( x2  +  y2  j 

Fifth-Order 

Aberration 

13 

3 

1 

3  x  - 1  2.y  (a-2  +  y2  j  + 1  O.y  ( A-2  +  y2  j 

Fifth-Order 

Aberration 

14 

3 

1 

3  y  - 1 2y  ( x2  +  y 2 )  + 1  Oy  ( x2  +  y2 ) 

Fifth-Order 

Aberration 

15 

3 

0 

-l  +  12(x2+y2)-30(x2+y2)2+20(x2+y2)3 

Fifth-Order 

Aberration 

16 

4 

4 

4  s-  2  7  4 

a*  -  6x  y  +  y 

Seventh-Order 

Aberration 

17 

4 

4 

4x3y  -  4xy3 

Seventh-Order 

Aberration 

18 

4 

3 

-4x3  +12at2  +5x3(x2  +y2)-15AT2(x2  +y2) 

Seventh-Order 

Aberration 

19 

4 

3 

-12x2y  +  4y3  +15x2y(x2  +y2)-5y3(x2  +y2) 

Seventh-Order 

Aberration 

Table  2.  First  21  Zernike  Polynomials  in  Cartesian  Coordinates  (From  Allen  2007) 


In  order  to  properly  fit  the  square  sensor  array  inside  of  a  unit  circle  the  SHWFS 
had  to  be  circumscribed  inside  of  the  unit  circle.  Assuming  an  NxN  square  sensor  grid, 
the  circle  was  created  with  a  radius  R  —  V2/V2  as  shown  in  Figure  8. 
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Figure  8.  SFIWFS  Grid  Circumscribed  inside  Unit  Disk  for  Zernike  Polynomial 

Evaluation 


Recalling  back  from  section  II.  A,  and  the  recent  explanation  of  Zemike  polynomials  it  is 
now  possible  to  define  any  wavefront  as  a  summation  of  m  modes,  multiplied  by  its 
corresponding  coefficient  a7 

m 

00-  y)  =  ^  ajZj(x,  y) 

7=1 

where  Zj  represents  the  jth  Zemike  mode.  Similarly  to  the  zonal  reconstruction  methods, 
it  is  possible  to  formulate  a  matrix  version  of  Zemike  mode  summation 


<p  —  [Z]  a 

—  *2 

where  0  is  N  length  vector,  a  is  a  m  length  vector  corresponding  to  the  mode 
coefficients,  and  [Z]  is  a  N“  by  m  size  matrix  corresponding  the  evaluation  of  each 
Zernike  mode  at  each  phase  location.  Using  the  pseudo  inverse  of  [Z]  the  coefficients,  a, 
can  be  solved  in  a  least  squares  solution. 

a  —  [Z]t0 

Recalling  that  Zernike  polynomials  are  orthonormal  over  the  entire  unit  disk, 
meaning  there  derivatives  exist  over  the  entire  unit  disk;  it  is  possible  to  fit  the  SHFWS 
slope  data  to  the  derivatives  of  the  Zemike  polynomials  (Southwell,  1980). 


dZjjx.y ) 
dx 
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5^ 


m 


I 


a; 


dZjjx.y ) 
dy 


The  same  matrix  formulation  can  apply  relating  the  SHWFS  slope  data  to  the 
Zernike  coefficients  through  a  least  squares  solution  as 

a  =  [dZ]fS 

where  dZ  is  a  matrix,  2N  by  m,  representing  the  derivatives  of  the  Zemike  polynomials 
for  both  x  and  y,  and  S  is  a  vector  of  length  2N2,  representing  the  x  and  y  slope 
measurements  of  the  SHWFS.  Once  the  coefficients  are  calculated  they  are  used  to 
directly  evaluate  the  phase  at  any  point  in  the  aperture,  or  corresponding  to  a  given 
lenslet.  The  higher  order  modes  are  related  to  higher  frequency  aberrations  in  a 
wavefront,  so  there  is  a  tradeoff  between  efficiency  and  accuracy,  when  deciding  how 
many  modes  to  use.  Modal  reconstruction  does  also  apply  a  form  of  smoothing  the  phase 
functions,  as  it  is  a  continuous  function,  which  can  introduce  some  errors  in  the  phase 
estimation. 


3.  Noise-Variance  Weighted  Complex  Exponential  Reconstructor 

It  was  discussed  previously  that  nulls  in  intensity  readings  on  the  SHWFS  lead  to 
the  presence  of  branch  points,  which  is  addressed  in  depth  by  Dr.  Fried  in  “Branch  Point 
problem  in  adaptive  optics.”  Branch  points  arise  as  an  issue  with  traditional  least  squares 
reconstructor,  like  zonal  or  modal  methods,  because  they  separate  the  phase  and 
amplitude  when  expressing  the  optical  field.  If  the  perturbations  are  expressed  as  a 
complex  phasor,  u  —  Aexp(i<p ),  composed  of  the  phase  and  amplitude,  (j)  and  A,  then 
branch  point  location  and  information  are  not  lost  during  the  reconstruction  (Fried,  2001). 
The  basic  concept  behind  the  NVWCER  lies  in  the  initial  formulation  of  how  to  express 
the  phases  and  phase  differences.  In  the  reconstruction  process  phases  are  represented  as 
phasors,  u  =  exp(i0),  and  phase  differences  are  represented  as  differential  phasors, 
A u  =  exp(iA0).  The  NVWCER  is  a  recursive,  multi-path  averaging  algorithm  which 
uses  multiple  paths,  from  an  initial  phasor  to  another  phasor,  to  best  estimate  the 
differences.  It  operates  on  a  Hudgin  based  geometry  grid  of  square  size  2N+1.  Like 
traditional  methods  multiple  paths  are  added  together  and  use  to  obtain  a  best  estimate  for 
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overall  difference,  however  with  phasors  one  multiplies  to  perform  the  equivalent  of 
adding  and  multiplies  by  the  complex  conjugate  for  subtraction. 

0i-02  ~  exp (i0i)  exp(-i02) 

0i  +  02  ~  exp(i01)  exp(i02) 


The  noise-variance  weighting  aspect  of  the  NVWCER  utilizes  the  relative  signal- 
to-noise-ratio  (SNR)  for  each  sensor  measurement  to  weight  the  corresponding  path 
versus  other  paths  for  the  calculation.  A  quick  review  of  statistics  leads  to  two  statements 
that  provide  the  basis  for  noise  variance  weighting  when  combining  several  noise  ladened 
quantities  (Fried,  2001).  When  noise  ladened  quantities  are  added  the  best  estimate  for 
the  value  of  the  sum  is  simply  the  sum  of  the  quantities,  and  the  noise  induced  variance  of 
the  sum  is  equal  to  sum  of  the  noise  induced  variances  of  each  of  the  quantity  (Fried, 
2001).  Secondly,  the  best  average  of  several  noise  ladened  quantities  is  the  sum  of  the 
ratios  of  the  quantity  divided  by  its  noise  induced  variance  -  that  sum  is  in  turn  divided 
by  the  sum  of  the  inverse  of  the  individual  noise  induced  variances  (Fried,  2001).  Using 
the  notation  of  for  noise  ladened  quantities  and  noise  induced  variance,  xn  and  a 
respectively: 


N 


N 


^ sum  ^  |  Xn  >  Osum  ^  '  °n 


%AVG  ~ 


n= 1 
yiV 

L>n=l  _2 
_ un 

y  N  _5_ 

£-m=  1  ^2 

Uri 


71  =  1 


°AVG  ~ 


y  N 

t-in=  1  rr-2 


Or, 


The  NVWCER  is  a  two  part  process,  first  the  differential  phasors  are  used  to 
reconstruct  the  wavefront  then  the  branch  points  are  analyzed  and  placed.  The 
reconstruction  is  a  three  step  process;  reduce,  solve,  and  rebuild  stages.  The  reduce  stage 
recursively  goes  from  a  2N  +  1  square  grid  to  a  2W_1  +  1  square  grid,  as  it  uses  several 
differential  phasors  and  calculates  a  single  noise-variance  weighted  differential  phasor 
corresponding  to  the  de-resolution,  until  a  2x2  grid  is  reached.  The  solve  stage  uses  the 
four  differential  phasors  from  the  last  step  of  the  reduce  stage,  along  with  an  assigned 
zero  phase  point,  to  calculate  the  phasor  at  each  corner.  The  rebuild  stage  uses  the 
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calculated  phasors  and  corresponding  differential  phasors  to  calculate  the  new  phasors 
working  from  a  2W  +  1  to  a  2N+1  +  1  square  grid  each  iteration,  until  the  original 
resolution  is  reach.  Each  stage  is  covered  in  depth  in  the  subsequent  sections 

a.  Reduce  Stage 

The  reduce  stage  of  the  algorithm  the  original  2N+1  square  grid  is  successively  reduced  to 
2n-1+1,  and  so  forth,  until  a  2  by  2  size  grid  is  obtained.  In  each  step  there  is  also  a 
corresponding  reduction  in  the  number  of  differential  phasors.  In  a  process  explained 
below,  the  differential  phasors,  and  accompanying  noise  induced  variances, 
corresponding  to  the  2N'1+1  square  grid  are  formed  from  a  noise-variance  weighted 
estimate  of  the  2N+1  grid  differential  phasors  and  variances  (Fried,  2001).  One  iteration 
of  the  reduce  stage  is  shown  in  Figure  9  where  grid  (a)  is  of  size  N=3,  and  is  reduced  to 
grid  (b)  of  size  N=2. 


Figure  9.  Fattice  Points  for  Reduce  Stage  of  NVWCER  (a)  N=3  (b)  N=2 

(From  Fried  2001) 


For  each  iteration  in  the  reduction  process  it  is  necessary  to  calculate  a 
value  for  the  differential  phasor  and  corresponding  noise  induced  variances  between 
adjacent  grid  point  on  the  2N1+1  grid.  The  differential  phasors,  for  either  x  or  y,  to  be 
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calculated  are  referred  to  as  xAu  and  yAu,  and  the  noise  induced  variances  are 
°xAU  and  ayAu-  F°r  notation  purposes  in  the  calculations  a,b,c,de,f,g,h,  and  i  refer  to 
differential  phasors  from  the  originating  lattice  in  each  step  of  the  reduction.  There  are 
six  possibilities  for  the  location  when  calculating  the  differential  phasor  for  the  reduced 
lattice.  The  first  four,  shown  in  Figure  10,  are  when  it  is  located  on  the  edge,  and  the 
second  two,  shown  in  Figure  11,  occur  on  the  interior  of  the  lattice.  The  lines  indicate 
differential  phasors  from  the  starting  lattice,  solid  circles  represent  grid  points  on  the 
reduced  lattice,  and  open  circles  represent  grid  points  on  starting  lattice. 


Figure  10.  Calculation  of  Differential  Phasors  for  Reduce  stage,  grid  points  located  on 

edge  of  lattice  (From  Fried  2001) 
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Figure  11.  Calculation  of  Differential  Phasors  for  Reduce  stage,  grid  points  located  on 

interior  of  lattice  (From  Fried  2001) 


When  the  differential  phasor  lies  between  two  grid  points  located  on  the 
edge  of  the  reduced  lattice,  the  following  four  pairs  of  equations  are  used.  It  is  important 
to  note  that  at  each  step  of  the  reduction  process,  the  corresponding  values  and  variances 
need  to  be  store  for  use  in  the  rebuild  stage.  For  Figures  10. a 
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For  Figure  10. b  the  following  equations  are  used 
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When  the  phasors  on  the  reduced  grid  lie  between  points  on  the  interior  of 
the  lattice  the  following  equations  are  used,  based  on  Figure  1 1.  For  Figure  1  l.a 
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This  provides  the  necessary  equations  for  carrying  out  the  reduce  stage  of 
NVWCER  from  a  2N+1  square  grid  to  a  2  by  2  square  grid.  Connecting  the  equations  with 
the  paths  laid  out  in  Figure  10  and  Figure  11,  along  with  the  mention  of  statistical 
averaging,  it  is  easy  to  follow  formulation  of  each  equation. 


b.  Solve  Stage 

The  reduce  stage  ends  with  four  differential  and  phasors  and  four  noise- 
induced  variances,  corresponding  to  the  outer  four  edges  of  the  lowest  resolution  lattice. 
If  a  phasor  value  of  0  is  assigned  to  one  of  the  corners  as  the  reference  point,  then  a  least 
squares  noise-variance  weighted  solution  for  the  phasor  at  each  corner  can  be  calculated. 
This  method  however  is  intended  for  use  with  phases  and  adapted  to  work  with  phasors, 
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so  a  method  is  implemented  to  work  with  phasors  directly  and  avoid  any  2n  ambiguities 
(Fried,  2013).  The  goal  is  to  determine  four  complex  value  phasors,  with  unity 
magnitude,  corresponding  to  the  four  corner  points,  of  the  grid  show  in  Figure  12. 


Figure  12.  Grid  Patter  for  Solve  Stage  of  NVWCER 


The  four  phasors,  analogous  to  A,  B,  C,  and  D  will  be  referred  to  as 
PA,PB,PC,  and  PDand  the  differential  phasors  from  the  reduce  stage  are 
uAba,uAcb,uAcd,  and  uAda.  The  relationship  between  the  phasors  and  measured 
differential  phasors  is  as  follows  (Fried,  2013): 


u^ba  — 


uAcb  ~  pT 


U^CD  ~  — 


U^DA  ~ 


Dr.  Fried  makes  note  though,  that  since  our  differential  phasors  have  a 
noise  induced  variance,  they  are  only  imperfectly  known,  so  a  path  for  example 
corresponding  from  A  to  D  yields 

u^da  ~  u^bau^cb(.u^cd') 

In  order  to  accurately  account  for  these  imperfect  measurements,  the  noise 


induced  variances  need  to  be  accounted  for  in  a  similar  fashion  to  the  statistical 
representation  in  the  reduce  stage.  First  an  absolute  phase  reference  is  created  by  assigned 
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a  value  for  position  A  of  P^  =  1.  One  example  path  calculation,  for  Point  C,  is  shown 
below  to  give  a  foundation  for  the  final  calculations.  Since  there  are  two  paths  from  A  to 
C,  it  is  required  to  compute  a  noise-variance  weighted  average  of  the  paths  (Fried,  2013). 
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making  note  that  the  variances  simply  add  along  the  path  to  get  the  corresponding  path 
noise  induced  variance  (oABC  =  oBA  +  °cb)-  Carrying  out  the  formulation  for  the 
remaining  paths  from  reference  point  A  to  the  remaining  points,  then  solving  in  tenns  of 
the  known  variances  and  differential  phasors  leads  to  the  final  equations  (Fried,  2013). 
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where  oALL  —  oBA  +  oBB  +  oBD  +  oBA  (Fried,  2013).  The  end  of  the  solve  stage  leaves 
he  algorithm  with  four  phasor  values  and  there  accompanying  noise  induce  variances. 
The  next  stage  uses  these  phasor  values  to  calculate  the  other  corresponding  grid  points. 


c.  Rebuild  Stage 

The  rebuild  stage  of  the  NVWCER  works  in  a  similar,  but  reverse  fashion, 
as  the  reduce  stage.  It  recursively  starts  with  a  2N+1  square  grid  of  phasor  values  and  uses 
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the  corresponding  differential  phasor  from  the  equivalent  step  of  the  reduce  stage  to 
calculate  the  a  noise-variance  weighted  value  for  the  phasors  on  a  2N+1+1  square  grid, 
constantly  repeating  until  the  original  lattice  size  is  achieved.  The  algorithm  first 
calculates  the  phasor  at  the  center  of  each  elemental  square,  then  the  middle  phasor  points 
on  the  edge  of  each  one  (Fried,  2013).  The  notations  A,  B,  C,  and  D  are  used  to  refer  to 
the  calculated  phasors  from  the  previous  step  of  the  rebuild  stage,  or  the  output  of  the 
solve  stage  if  it  is  the  first  iteration.  The  use  of  x#  and  y#  are  used  to  represent  the 
differential  phasors  from  the  appropriate  step  of  the  reduce  stage,  with  the  same  amount 
of  current  grid  points.  The  equation  for  the  phasor,  u,  and  its  noise  induced  variance, 
at  the  center  of  the  elemental  square  follow  the  notation  shown  in  Figure  13,  where  the 
solid  circle  inside  an  open  circle  represents  the  phasor  to  be  calculated 


o  i-o  . ;•  o  .  i-  o  . i-O 


•  o  . o  .  o  . i-  o  . i-  o  ••• 

Figure  13.  Calculation  for  center  phasor  of  each  elemental  square  in  Rebuild  stage  (From 

Fried  2001). 
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Once  the  phasors  at  the  center  of  all  the  unit  squares  are  calculated,  the 
next  step  is  to  calculate  the  phasors  at  the  middle  of  the  edges  of  the  unit  squares.  When 
the  edge  of  the  unit  square  lies  on  the  interior  of  the  lattice  then  the  following  equations 
are  used,  representing  the  notation  in  Figure  14.  The  sold  black  circles  represent  phasors 
from  the  previous  steps  of  the  rebuild  stage,  the  small  dots  with  open  circles  around  them 
represent  the  calculated  phasors  at  the  center  of  each  elemental  square,  and  the  small 
black  circles  are  the  phasors  to  be  calculated. 
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Figure  14.  Calculation  for  Phasors  at  the  middle  of  edges  of  elementary  squares,  located 

on  interior  of  lattice  (From  Fried  2001) 


The  calculations  apply  to  both  setups  in  Figure  14,  either  edge  on  top  or 
bottom  of  elemental  square  (left)  or  on  the  left  or  right  of  elemental  square  (right). 
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If  the  phasor  is  located  on  the  edge  of  the  lattice  then  there  are  four  possibilities; 
either  on  the  sides  of  the  lattice  or  on  the  top  or  bottom  of  the  lattice.  Figure  15  shows  the 
notation  for  all  four  possibilities,  and  the  same  symbols  apply  from  Figure  14. 
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(a)  (b) 


Figure  15.  Phasor  calculation  setup  for  Rebuild  stage;  phasor  located  on  edge  of  lattice 

(From  Fried  2001) 


Dependent  on  where  the  phasor  is  located  differing  set  of  equations  apply. 

For  a  phasor  located  on  the  bottom  of  the  lattice,  show  in  Figure  15. a, 
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For  phasors  located  on  the  top  as  show  in  Figure  15.b; 
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For  phasors  located  on  the  left  edge  as  shown  in  Figure  15.c; 
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For  phasors  located  on  the  right  edge  as  show  in  Figure  15.d; 
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With  the  three  stages  presented,  along  with  the  supporting  equations  and 
notation,  it  is  possible  to  carry  out  the  full  extent  of  the  NVWCER.  However,  this 
algorithm  is  only  partially  complete,  because  it  has  done  nothing  for  unwrapping  or 
branch  point  placement.  This  is  purely  a  calculation  of  the  phase  in  the  complex  domain, 
with  no  consideration  of  branch  points  yet.  A  second  part  of  the  algorithm  will  take  the 
complex  phasors  attempt  to  identify  and  place  the  branch  points,  as  well  as  provide  an 
unwrapped  and  continuous  phase  in  the  real  domain. 


d.  Phase  Unwrapping  and  Branch  Point  Analysis 

Give  the  output  of  phasors  from  the  rebuild  stage  of  the  NVWCER,  there 
are  two  main  issues.  First  it  is  in  the  complex  domain  as  phasors  and  not  phases,  which 
are  required  for  an  AO  system.  Secondly,  the  system  has  not  accounted  for  branch  points, 
in  which  there  will  be  large  2%  discontinuities.  This  “smoothing”  algorithm,  presented  by 
Dr.  Fried,  takes  the  phasor  data  and  extracts  a  continuous,  unwrapped  phase  function  that 
has  ideally  placed  branch  points  close  together  and  in  areas  of  low  intensity  as  to 
minimize  the  effects  on  the  AO  system. 
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The  algorithm  first  starts  by  calculating  the  differential  phasors  from  the 
output  phasors  of  the  rebuild  stage;  these  ideally  are  the  same  measurements  as  the  initial 
stage,  however  with  noise  considerations  taken  (Fried,  2001).  These  differential  phasors 
are  summed  up  around  in  each  unit  circle  to  identify  any  branch  points.  The  next  step 
uses  the  location  of  the  branch  points  to  calculate  the  hidden  phase,  and  remove  that  from 
the  phasors,  ideally  leaving  a  scalar,  branch  point  free  phase  function.  Using  the  scalar 
phase  function,  which  ideally  has  no  branch  points,  the  phase  can  be  calculated  by  simply 
setting  a  reference  zero  phase  point  and  just  adding  the  differential  phasors  to  calculate 
each  phasor.  Finally,  the  total  phase  function  is  created  as  a  sum  of  the  hidden  phase  and 
scalar  phase  function,  which  has  accounted  for  branch  points  and  placed  them  preferably 
in  locations  of  low  intensity. 

First,  the  phase  differences  are  calculated  across  the  entire  lattice 

according  to 

xA(j)ij  =  arg(ui+1ju*j ) 

yA<pi,j  =  arg(uUjulJ+ 1) 

and  then  the  curl,  C(i.j),  is  calculated  about  each  unit  square,  or  sub-aperture,  according 
to 

C(ij')  =  xAtptj  +  yA(pi+lj  -  xAc/)iJ+1  -  yA<pUj 
Once  the  curls  are  calculated  then  the  branch  points  can  be  located.  Assuming  the  range 
of  the  phases  is  from  -71  to  +71,  then  that  is  used  for  the  reference  of  finding  branch  points. 
A  positive  branch  point  is  denoted  by  P  =  {(ip,jp)}  and  negative  branch  points  by 
N  —  {(inUn))  (Fried,  2001).  The  criteria  for  being  classified  a  branch  points  is  as 
follows: 

( J-p’jp )  £  P  if  Cif-pijp')  -'>  d"7r 
( ip.jp )  £  P  if  C{ip,jp)  <  +7T 
C in’jn)  ^  P  if  P(jn>jn)  <~- 
iirvjn )  £  P  if  C(in,jn )  ^  ~Tt 

Since  the  indexing  for  each  curl  function,  C(i,j),  is  indexed  off  of  the 
bottom  left  corner  of  each  unit  square,  then  it  is  necessary  to  account  for  the  location  of 
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the  branch  point  being  located  at  the  center.  The  actual  locations  of  the  branch  points  are 
denoted  pp  =  (f, ip,r]jp )  andpn  =  ($in,Vjn),  where 

Zb  =\{xiP  +  %«)'  vJp  =  \{y]p  +  y,P+1) 

^  hi  =  2  ^X'n  +  Xhi+i)’  Vjn  =  2  O^n  +  ^/n+l)  ' 

Once  the  curls  around  each  unit  square  are  calculated,  and  all  of  the 
branch  points  are  detected,  the  total  amount  of  branch  points  is  tallied.  If  there  are  more 
positive  (negative)  branch  points  then  negative  (positive),  an  artificial  negative  (positive) 
branch  point  is  placed  at  the  maximum  (minimum)  of  the  potential  field,  V(i,j).  This 
process  is  repeated  until  the  number  of  positive  and  negative  branch  points  are  equal. 
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Then  the  hidden  phase  is  calculated  and  removed  from  the  original  phasor 
lattice  to  create  a  scalar  function,  uL  j,  free  of  branch  points. 
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The  scalar  function  is  free  of  branch  points,  so  the  phase  can  be  calculated  through  a 
simple  summation  through  the  paths  of  differentials.  When  adding  up  the  scalar  phase 
differences,  a  reference  value  of  0  is  set  at  an  arbitrary  point. 

xHi,j  =  argfa+ijulj ) 

=  arflf(uu+iu-j) 
t'=i-i  j'=j- 1 

tpscalari-i’ f)  ~  1  xAfoj  +  I  yAcptji 

i'= 1  j'= 1 

With  the  scalar  phase  and  hidden  phase  calculated,  the  total  phase  is  calculated  as  a 
summation  of  the  hidden  and  scalar  phase  functions. 


0(hf)  0hid(hf)  T  (pscalari-i’ ]) 

This  process  creates  a  phase  which,  as  best  as  possible,  is  unwrapped  with 
branch  points  placed  is  areas  of  lowest  optical  intensity.  It  is  important  to  note  again  that 


40 


this  process  is  defined  for  Hudgin  geometry.  The  next  section  discusses  how  to  properly 
adapt  this  algorithm  for  a  fried-based  geometry,  making  it  compatible  with  a  SHWFS. 

e.  Adaptation  to  Fried  Geometry 

Previously  mentioned,  the  SHFWS  operates  on  either  Fried  or  Southwell 
geometry,  but  the  NVWCER  operates  on  Hudgin  geometry.  The  simplest  way  is  to 
convert  the  SHFWS  data  into  two  interleaved  Hudgin  lattices,  carry  out  the  NVWCER 
process  on  each  of  them,  and  then  interleave  them  back  together.  Looking  at  a  basic 
Hudgin  lattice,  as  shown  in  Figure  16,  the  black  dots  represent  the  points  on  one  Hudgin 
lattice,  and  the  white  dots,  the  other  lattice. 


Figure  16.  Fried-to-Hudgin  Geometry  Conversion  (From  Fried  2001) 

The  phase  differences  in  the  corresponding  Hudgin  lattices  show  in 
Figure  17  have  phase  differences  which  are  the  difference  of  diagonally  opposite  phases 
in  a  given  supaperture.  In  Fried  geometry,  the  x  or  y  slope  measurements  of  a  SHWFS 
correspond  to  the  difference  of  the  average  of  the  two  phases  on  each  side,  for  x,  or  top 
and  bottom  for  y,  divided  by  the  distance  (Fried,  2001).  Mathematically,  the  diagonal 
phase  differences,  as  shown  in  the  Hudgin  lattices,  correspond  to  the  sum  or  difference  of 
those  same  x  and  y  slope  measurements  from  the  SHWFS  for  that  given  supaperature 
(Fried,  2001). 
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Figure  17.  Corresponding  Hudgin  Geometries  for  a  SHFWS  Fried  Geometry  (From  Fried 

2001) 

Given  both  of  the  lattices  in  Figure  17,  values  not  corresponding  to  the 
original  lattice  are  given  a  phasor  value  of  0  and  a  variance  value  of  infinity,  this  in 
essence  allows  the  NVWCER  to  operate  without  actually  using  them  in  subsequent 
calculations.  Once  the  NVWCER  is  ran  on  both  lattices  the  phase  points  are  then  inserted 
back  into  the  original  positions  to  form  a  Fried-geometry  lattice.  Since  there  is  no 
absolute  reference  phase  point,  only  a  relative  reference  phase  for  each  individual  lattice, 


there  is  a  chance  they  are  out  of  phase  with  each  other.  In  order  to  bring  them  back  in 
phase,  the  phases  corresponding  to  the  white  dots  in  Figure  16  are  compared  with  an 
average  of  the  four  adjacent  grid  points  to  calculate  the  average  phase  shift.  This  phase 
shift  is  then  applied  to  the  lattice. 


First  the  noise-variance  weighted  average  of  the  four  adjacent  grid  points 
to  each  white  grid  point  from  Figure  16  is  calculated.  The  phases  of  the  un-synced  lattice 
are  referred  to  as  u(i,j ). 


where 


u(i  -  1 u(i  +  1 ,  u(i,j  -  1)  u(i,j  +  1) 

■\  i  9  r  ■  .  *\  ”r  o  r  •  •  ^  \  T 


fi(ij)  = 


o2(i-l,j)  o2  (i  +  1,  j)  o2(i,j-l)  <J2(i,  j  +  1) 


for  i  +  j  od 


T\  + 


T\  + 


+ 


o2(i-l,j)  o2(i  +  1,  j)  <J2(i,  j  —  1)  a2(i,j  +  l) 


Then  the  weighted  sum  of  the  products  of  the  averaged  value,  and 

the  conjugate  of  the  original  un-synced  phases,  u(i,j ),  is  found  to  calculate  the  weighted 
average  phase  shift  (Fried,  2001). 


0(iJ)  =  arg{  ^ 


...  S2{i,j) 

i+j  odd 

a2(i,;)  =  52(i,j )  +  a2(i,j )  ,  for  i  +  j  odd 


Once  the  average  phase  shift  is  applied  the  final  phasor  values,  u(i.j),  are 
calculated  and  the  phase  shift  is  applied  to  correct  Hudgin  lattice  points. 


u{i.j) 


f  u(i,j )  for  i  +  j  odd 

(exp (i0)  u(i.j )  for  i  +  j  even 


This  final  phasor  lattice  can  then  be  applied  to  the  branch  point  analysis 
and  unwrapping  algorithm.  It  is  important  to  note  that  due  to  the  interleaving  of  two 
lattices,  and  a  finite  estimate  of  the  phase  shift,  there  is  subsequent  error  that  can 
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propagate  forward  as  a  result.  The  attempt  to  re-sync  the  phases  can  lead  to  minor  “egg- 
crating”  effects  due  to  areas  of  minor  phase  shift  discrepancies  (Fried,  2001) 
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V.  ANALYSIS 


A  MATLAB  simulation  environment  is  created,  to  analyze  the  different 
reconstruction  methods  with  varying  Rytov  variances,  SNR,  and  sensor  grid  size.  It  is 
important  to  note  when  random  noise  is  added  to  generate  the  desired  SNR,  the 
corresponding  SHFWS  output  is  simultaneously  sent  to  all  three  reconstructors.  This 
assures  comparison  is  made  using  the  exact  same  sensor  data.  The  principal  measurement 
for  the  quality  of  reconstruction  is  known  as  the  Strehl  ratio.  It  is  a  weighted 
measurement  of  the  difference  between  original  phase  and  reconstructed  phase  at  each 
point.  The  Strehl  ratio  is  described  by  the  following  equation, 


5  = 


1  [li= i  uorigianl(i,j)u*rec(i,j )  ]  |' 

[Y.i=l'L'j=l\Uoriginal(Uj)\]2 


where  it, 


origia.nl 


is  the  complex  function  of  the  original  phase  and  amplitude  sampled  at 


the  reconstructed  points  and  urec  is  the  complex  function  of  the  reconstructed  phase  with 
unity  amplitude. 


uoriginal(i>j )  A(i,j~)e  ^ 

UreciUj )  =  e^recVJ) 

Looking  at  the  equation  for  the  Strehl  ratio,  where  the  value  of  the  original  is 
multiplied  by  the  complex  conjugate  of  the  reconstructed  that,  it  is  essentially  calculating 
the  difference  in  phase.  The  closer  the  values  are,  the  closer  the  values  of  the  exponential 
are  driven  to  zero,  leading  to  a  value  of  one. 


All  reconstructors  are  evaluated  using  all  possible  combinations  of  the  following 
factors:  SNR  (10,  60,  and  200),  Number  of  sensors  along  one  side  (16,  32,  and  64),  and 
Rytov  variances  (0.109,  0.349,  and  0.567).  In  each  section  of  the  analysis  only  certain 
results  are  discussed,  in  order  to  obtain  an  accurate  analysis.  A  table  of  results  for  all 
combinations  is  included  as  well 


A.  SOURCE  PHASE  AND  INTENSITY  DATA 

The  data  used  in  testing  and  analysis  of  the  various  reconstruction  algorithms  was 
generated  using  a  program  known  as  WaveProp.  This  is  a  proprietary  MATLAB  program 
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which  simulates  the  propagation  of  a  beam  through  a  turbulent  atmosphere,  with  given 
parameters  like  distance,  turbulence,  wavelength,  and  size.  The  data  used  in  this 
simulation  had  the  parameters  of  a  propagation  distance  of  4km,  a  wavelength  of  one 
micron  (lpm),  and  a  sensor  side  length  of  0.3m.  The  input  for  WaveProp  is  ,  but  given 
the  path  length  and  wavelength,  the  equation  for  Rytov  variance  can  be  used  to  calculate 
the  proper  C%  for  a  desired  Rytov  number. 

aj  =  1.23Cn2  (y)^ 

The  following  are  images  of  the  phase  and  amplitude  data. 


Amplitude  x  10*3  Phase 


Figure  18.  Simulation  data  for  Rytov  Variance  of  0. 109 
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Phase 


Figure  19.  Simulation  data  for  Rytov  Variance  of  0.349 


Phase 
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Figure  20.  Simulation  data  for  Rytov  Variance  of  0.567 

The  phase  data  presented  is  wrapped,  meaning  it  is  on  the  domain  of  zero  to  2ti. 
The  interest  of  reconstruction  is  generating  an  unwrapped  phase,  meaning  it  exist  on  the 
domain  of  —  oo  to  +  oo,  so  that  it  is  as  close  to  a  continuous  as  possible.  This  relates  to  the 
need  for  generating  commands  for  the  deformable  mirror,  which  is  continuous  as  well 
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and  cannot  create  step  discontinuities.  When  results  are  presented  both  the  wrapped  and 
unwrapped  phase  will  be  shown,  however  in  the  Strehl  ratio  calculation  it  has  no  effect 
due  to  the  domain  a  complex  exponential  repeating  every  2n. 

B.  INITIAL  TEST 

1.  Variation  in  Number  of  Sensors 

The  number  of  sensors  is  evaluated  at  16,  32,  and  64  along  one  for  each 
reconstructor.  The  data  presented  in  this  section  for  analysis  has  an  SNR  of  60,  and  is 
analyzed  for  a  Rytov  variance  0.567. The  primary  concern  is  to  see  if  there  is  a  limit  to  the 
number  of  sensors  needed  for  a  relatively  high  Rytov  Variance. 
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Modal  Wrapped  Strehl=0.067748  SNR=60N=16 


Figure  22.  Modal  Reconstruction.  S=0.068,  SNR=60,  and  N=16 


Zonal  Wrapped  Strehl=0. 12596  SNR=60N=16 


Figure  23.  Zonal  Reconstruction.  S=0.126,  SNR=60,  and  N=16 
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The  reconstruction  of  a  high  Rytov  variance  of  0.567  is  difficult  with  a  small 
amount  of  sensors.  The  Modal  reconstruction  performs  the  worst  with  a  Strehl  of  0.068, 
which  is  a  combination  of  the  smoothing  that  is  inherent  to  the  modal  reconstruction 
process  and  a  low  sensor  count.  The  Zonal  and  NVWCER  performed  nearly  identical, 
with  StrehTs  of  0.126  and  0.146,  respectively.  Zonal  reconstruction  is  a  method  proven  to 
work  well  with  lower  disturbance  regimes,  and  given  a  lower  sensor  count  the  NVWCER 
has  some  difficulties  properly  identify  the  branch  point  locations.  The  following  set  of 
data  is  for  32  sensors  along  a  side. 


CER  Wrapped  Strehl=0.31224  SNR=60N=32 


5  10  15  20  25  30 


Figure  24.  NVWCER  Reconstruction.  S=0.3 12,  SNR=60,  N=32 
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Modal  Wrapped  Strehl=0. 1 1 428  SNR=60N=32 


Figure  25.  Modal  reconstruction.  S=0. 1 14,  SNR=60,  and  N=32 


Zonal  Wrapped  Strehl=0.23784  SNR=60N=32 


5  10  15  20  25  30 

Figure  26.  Zonal  Reconstruction.  S=0.238,  SNR=60,  and  N=32 
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When  the  number  of  sensors  is  increased  by  a  factor  of  2,  to  32,  the  Strehl 
increases  for  all  three  reconstruction  methods  as  anticipated.  The  Modal  again  perfonns 
the  worst  of  all  three  reconstructors.  The  NVWCER  outperforms  the  Zonal  reconstructor 
by  a  larger  margin  this  time,  with  StrehTs  of  0.312  and  0.238,  respectively.  Now  that  the 
NVWCER  is  receiving  adequate  data,  it  is  able  to  perfonn  better  in  the  presence  of  high 
turbulence,  as  designed,  then  other  methods.  It  is  also  important  to  note  that  in  Figure  24, 
the  checkerboard  pattern  is  beginning  to  emerge,  which  was  attributed  to  the  process  of 
the  NVWCER  for  a  SHFWS  design  needing  to  bring  the  two  lattices  into  phase  with  each 
other.  The  last  set  of  data  is  for  a  grid  of  64  sensors  along  a  side. 


CER  Wrapped  Strehl=0.29749  SNR=60N=64 


10  20  30  40  50  60 

Figure  27.  NVWCER  Reconstruction.  S=0.297,  SNR=60,  and  N=64 
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Modal  Wrapped  Strehl=0. 1 1 573  SNR=60N=64 


10  20  30  40  50  60 

Figure  28.  Modal  Reconstruction.  S=0.1 16,  SNR=60,  and  N=64 


Zonal  Wrapped  Strehl=0. 20619  SNR=60N=64 


10  20  30  40  50  60 


Figure  29.  Zonal  Reconstruction.  S=0.206,  SNR=60,  and  N=64 


The  reconstructors  all  perforin  in  the  same  manner  as  the  previous  two  sections 


with  modal  performing  the  worst  and  NVWCER  perfonned  the  best.  However,  there  was 
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a  decrease  in  the  Strehl  ratio  of  both  the  NVWCER  and  zonal  reconstructor  of  0.2  and 
0.3,  respectively,  which  would  lead  one  to  believe  that  it  is  possible  to  over  sample  the 
data.  Possible  outlying  factors  could  also  be  the  location  density  of  the  random  noise 
added  to  the  system.  While  the  SHWFS  data  is  sent  to  all  three  reconstructors 
simultaneously  for  any  trial,  if  any  of  the  factors  are  changed  a  new  random  matrix  of 
noise  data  is  generated;  the  SNR  is  on  a  measurement  defined  over  the  entirety  of  the 
image. 

The  variation  in  the  number  of  sensors  has  shown  a  trend  that  increasing  the 
number  of  sensors  to  a  threshold  will  have  a  positive  correlation  on  the  reconstruction 
process  and  the  Strehl  ratio.  The  NVWCER  has  shown  it  requires  a  minimum  amount  of 
sensors  to  outperfonn  the  zonal  reconstructor  to  a  significant  degree,  which  is  attributed 
to  the  fact  it  needs  enough  information  in  order  to  locate  and  handle  the  branch  points. 
The  next  section  with  analyze  the  variation  of  SNR. 

2.  Variation  in  SNR 

The  SNR  is  used  to  represent  the  relative  capabilities  of  the  sensor  its  self. 
Depending  on  the  SHWFS  sensitivity  the  SNR  can  vary  to  large  degree.  An  SNR  of  10  is 
considered  relatively  low,  while  a  value  of  60  is  more  realistic,  and  200  is  a  sensors 
susceptible  to  very  little  noise.  The  reconstructors  will  be  tested  with  all  three  SNR’s.  A 
Rytov  variance  of  0.349  will  be  used  and  the  number  of  sensors  will  be  set  to  32  on  a 
side. 
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CER  Wrapped  Strehl=0.14412  SNR=10N=32 


5  10  15  20  25  30 


Figure  30.  NVWCER  Reconstructor.  S=0. 144,  SNR=10,  N=32 


Modal  Wrapped  Strehl=0. 16957  SNR=10N=32 


5  10  15  20  25  30 


Figure  31.  Modal  Reconstruction.  S=0.169,  SNR=10,  and  N=32 
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Zonal  Wrapped  Strehl=0.2956  SNR=10N=32 


5  10  15  20  25  30 


Figure  32.  Zonal  Reconstruction.  S=0.296,  SNR=10,  and  N=32 

With  such  a  low  SNR  all  of  the  reconstructors  obtain  a  relatively  low  Strehl  ratio. 
The  zonal  reconstructor  performs  the  best,  with  a  value  of  0.296,  the  modal  process  with 
a  Strehl  of  0.169,  and  NVWCER  with  a  Strehl  of  0.144.  The  modal  reconstructors 
inherent  smoothing  assisted  with  the  high  amount  of  noise.  The  NVWCER  is  designed  to 
handle  noisy  measurements  however,  with  an  SNR  of  10,  all  paths  of  a  high  amount  of 
noise.  The  noise  may  also  be  masking  much  of  the  detailed  information  needed  to  located 
and  handle  branch  points.  The  SNR  is  now  increased  to  60. 
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CER  Wrapped  Strehl=0. 73742  SNR=60N=32 


5  10  15  20  25  30 


Figure  33.  NVWER.  S=0.737,  SNR=60,  and  N=32 


Modal  Wrapped  Strehl=0.3558  SNR=60N=32 
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Figure  34.  Modal  Reconstruction.  S=0.356,  SNR=60,  and  N=32 
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Zonal  Wrapped  Strehl=0.54397  SNR=60N=32 


5  10  15  20  25  30 

Figure  35.  Zonal  Reconstruction.  S=0.544,  SNR=60,  and  N=32 


When  the  SNR  is  increased  to  60,  representing  a  relatively  sensitive  SHWFS  but 
still  subjected  to  noise,  all  of  the  reconstructors  increase  in  Strehl  ratio.  The  NVWCER 
performed  the  best  with  a  Strehl  of  0.737,  with  the  zonal  performing  second  with  a  Strehl 
of  0.544,  and  modal  performing  the  worst  with  a  Strehl  of  0.356.  With  the  NVWCER 
receiving  less  noisy  data,  it  is  able  to  perform  as  expected  in  the  presence  of  deep 
turbulence.  The  last  section  is  an  SNR  of  200,  which  represents  a  nearly  ideal  SHWFS. 
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CER  Wrapped  Strehl=0.76316  SNR=200N=32 


Modal  Wrapped  Strehl=0.35381  SNR=200N=32 


5  10  15  20  25  30 


Figure  37.  Modal  Reconstruction.  S=0.354,  SNR=200,  and  N=32 
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Zonal  Wrapped  Strehl=0.53634  SNR=200N=32 


5  10  15  20  25  30 

Figure  38.  Zonal  Reconstruction.  S=0.536,  SNR=200,  and  N=32 


A  SNR  of  200  yielded  minimal  differences  from  a  value  of  60,  with  only  a  slight 
increase  in  the  Strehl  for  the  NVWCER.  The  reconstructors  already  performed  decently 
with  an  SNR  of  60,  and  the  noise  was  not  the  driving  factor  at  that  point.  It  was  just  the 
ability  of  the  sensors  to  handle  the  SHWFS  data  in  the  presence  of  branch  points.  Since 
the  NVWCER  is  designed  to  handle  branch  points,  it  makes  sense  that  it  would  be  the 
one  to  gain  the  most  from  a  higher  SNR,  then  the  other  reconstructors. 

The  variations  in  both  SNR  and  number  of  sensors  have  yielded  interesting 
results,  and  from  them,  the  data  presented  to  compare  all  of  the  reconstructors  will  use  an 
SNR  of  60  and  an  array  of  sensors  with  length  32  on  a  side. 

C.  EVALUATION  OF  RECONSTRUCTORS 

The  reconstructors  are  evaluated  for  each  Rytov  variance,  with  a  constant  SNR 
and  number  of  sensors.  The  zonal  and  modal  reconstructors  are  designed  to  work  in 
environments  with  little  to  no  turbulence,  while  the  NVWCER  is  designed  specifically  to 
handle  branch  points,  as  a  result  of  high  turbulence.  The  0.109  Rytov  variance  represents 
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an  area  of  low  turbulence,  so  all  reconstructors  should  perform  similarly,  while  the 
variances  of  0.349  and  0.567  are  areas  of  higher  turbulence,  where  the  NVWCER  should 
outperform  the  others. 

1.  Rytov  Variance  of  0.109 
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Zonal  Wrapped  Strehl=0.98835  SNR=60N=32 
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Figure  39.  Reconstructor  Comparison  for  Rytov  of  0. 109.  SNR=60  and  N=32 


With  a  low  about  of  turbulence,  the  reconstructors  perfonned  as  predicted,  with 
the  zonal  reconstructors  and  NVWCER  performing  nearly  the  same  with  an  almost  a 
perfect  Strehl  ratio;  0.977  and  0.988,  respectively.  The  modal  reconstructor  is  still  a 
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victim  of  its  inherent  smoothing  leading  to  a  lower,  0.862,  Strehl  ratio,  with  a  limited 
number  of  modes  used. 

2.  Rytov  Variance  of  0.349 

CER  Wrapped  Strehl=0.73742  SNR=60N=32  Modal  Wrapped  Strehl=0.3558  SNR=60N=32 
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Figure  40.  Reconstructor  Comparison  for  Rytov  of  0.349.  SNR=60  and  N=32 

When  the  turbulence  is  increased  the  effects  of  each  reconstructor  are  more 
prominent.  The  modal  performs  significantly  less,  with  a  Strehl  of  0.356.  The  NVWCER 
is  outperforms  the  zonal  reconstructor,  with  Strehl  ratios  of  0.737  and  0.544,  respectively, 
due  the  greater  presence  of  branch  points  and  branch  cuts.  Since  the  zonal  reconstructor 
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does  not  account  for  branch  points,  the  more  of  them  present,  the  greater  the  error  of  the 
reconstructor. 


3.  Rytov  Variance  of  0.567 

CER  Wrapped  Strehl=0.31224  SNR=60N=32 


10  20  30 


Modal  Wrapped  Strehl=0. 11428  SNR=60N=32 


Zonal  Wrapped  Strehl=0.23784  SNR=60N=32 
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Figure  41.  Reconstructor  Comparison  for  Rytov  of  0.567.  SNR=60  and  N=32 


When  the  turbulence  was  increased  again,  to  a  Rytov  variance  of  0.567,  all  three 
reconstructors  drop  significantly  in  the  Strehl  ratio.  The  modal  reconstructor  perfonned 
the  worst  again,  with  a  Strehl  of  0.114,  while  the  NVWCER  outperfonned  the  zonal 
reconstructor  with  Strehl’s  of  0.312  and  0.238  accordingly.  The  ability  of  any 
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reconstructor  greatly  falls  off  when  the  turbulence  increases  beyond  a  certain  threshold, 
and  it  is  no  longer  an  issue  of  optical  correction,  rather  an  inappropriate  operational 
environment. 


D.  FINAL  RESULTS 

Table  3  is  the  final  results  of  all  simulations  ran  for  complete  comparison  and 
evaluation  of  all  three  reconstructors.  All  permutations  of  possible  SNR  (10,60  and  200), 
Rytov  variances  (0.109,  0.349,  and  0.567),  and  number  of  sensors  per  side  (16,  32,  and 
64)  were  ran  and  evaluated. 


Strehl  Ratio 

Sensors 

16 

32 

64 

SNR 

10 

60 

200 

10 

60 

200 

10 

60 

200 

Zonal 

0.109 

0.8461 

0.9699 

0.9649 

0.8645 

0.9884 

0.9897 

0.8350 

0.9887 

0.9968 

0.349 

0.1812 

0.4055 

0.4594 

0.2304 

.5292 

0.5403 

0.2205 

0.7015 

0.6178 

0.567 

0.1201 

0.1405 

0.1330 

0.1515 

0.2035 

0.2052 

0.0920 

0.2295 

0.2312 

Modal 

0.109 

0.7800 

0.8597 

0.8574 

0.7796 

0.8624 

0.8626 

0.7863 

0.8560 

0.8648 

0.349 

0.1273 

0.2886 

0.3277 

0.1561 

0.3445 

0.3513 

0.1402 

0.4370 

0.3929 

0.567 

0.0478 

0.0809 

0.0714 

0.0943 

0.1057 

0.1082 

0.0437 

0.1237 

0.1208 

NVWCER 

0.109 

0.6016 

0.9538 

0.9584 

0.6331 

0.9796 

0.9873 

0.2539 

0.9805 

0.9915 

0.349 

0.0187 

0.5557 

0.5665 

0.2063 

0.7831 

0.7724 

0.1328 

0.8737 

0.8907 

0.567 

0.0877 

0.1609 

0.1882 

0.1093 

0.3696 

0.4017 

0.2093 

0.6266 

0.3942 

Rytov 

Variance 

Table  3.  Results  for  all  Reconstruction  Simulations 
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VI.  CONCLUSION 


This  research  was  aimed  at  computationally  evaluating  the  perfonnance  of  two 
well-known  weave  front  reconstruction  methods,  zonal  and  modal,  against  a  modern 
reconstruction  method,  a  Noise- Variance  Weighted  Complex  Exponential  Reconstructor, 
created  by  Dr.  David  Fried,  in  the  presence  of  deep  turbulence.  The  NVWCER  is 
designed  to  reconstruct  using  a  complex  exponential,  so  it  can  utilize  both  the  hidden  and 
scalar  phase,  allowing  it  better  handle  branch  points.  The  reconstructors  were  evaluates 
using  simulated  SHFWS  data  using  various  combinations  of  Rytov  variances,  SNR,  and 
number  of  sensors. 

After  testing  an  evaluation  the  NVWCER  has  shown  its  ability  to  outperform  the 
modal  and  zonal  reconstruction  methods  in  the  presence  of  deep  turbulence.  At  a  Rytov 
variance  of  0.109,  representing  weak  turbulence,  all  three  reconstructors  perfonned  well. 
When  the  turbulence  level  was  increased  the  NVWCER  continuously  perfonned  better 
than  both  other  reconstructors.  The  zonal  reconstruction  method  typically  outperformed 
the  modal  as  well,  due  to  its  inherent  smoothing.  There  are  some  outlying  effects  such  as 
when  the  number  of  sensors  was  increased  greatly;  it  sometimes  had  a  negative  effect  on 
the  NVWCER.  The  reasoning  for  this  can  only  be  postulated  presently,  as  the  sample 
points  being  smaller  than  the  branch  point  themselves.  The  effects  of  “egg-crating”  also 
were  a  source  of  error  for  the  NVWCER,  which  was  due  the  use  of  two  lattices  and 
bringing  them  into  phase  with  each  other,  in  order  to  operate  with  the  SHWFS  geometry. 

This  research  shows  promise  for  the  use  of  a  NVECER  algorithm  in  environments 
with  deep  turbulence,  like  a  maritime  environment.  Two  areas  of  concern  and  interest  for 
future  work  involve  the  computational  efficiency  of  the  NVWCER  as  well  as  a  different 
sensor  pairing.  Currently  the  NVWCER  is  an  iterative  process,  using  many  for  loops  and 
minimal  matrix  operations,  and  it  would  be  useful  to  have  some  form  of  a  parallel  process 
designed.  The  most  important  issue  with  the  NVWCER,  and  its  use  with  a  SHFWS,  is  the 
differences  in  geometries.  The  NVWCER  is  designed  for  a  Hudgin  geometry,  and  later 
altered  to  work  with  Fried  or  Southwell  geometry,  however  that  introduces  the  lattice 

phase  issues,  which  in  turn  induce  error.  Either  a  shearing  interferometer  needs  to  be 
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paired  with  the  NVWCER,  or  the  algorithm  needs  to  be  adapted  to  work  directly  with 
Fried  or  Southwell  geometry. 
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APPENDIX 


A.  INITIALIZATION  AND  SIMULATION  CODE 

close  all; 
clear  all; 
clc; 

%%  Phase  Loading 

%  Choose  one  for  random  phase  generation  or  2  for  MyPhs  data 
vv  =  3  ; 

if  vv==3 

load ('amp  n  phase  0.567  l.mat'); 
phase2pi=phase2pi+pi ; 

N=9; 

phil=phase2pi ; 

end 

SNR0=60 ; 

NN=2 A6; 

NoiseFlag=l;  %1  adds  noise.  0  does  NOT  add  noise 
% [Fx, Fy, Magnitudes, SNR]  = 

slope  WgtAvg (phase2pi ,  ones (size (phase2pi) )  ,  NN,  SNRO ,  NoiseFlag) 
[Fx, Fy, Magnitudes , SNR]  = 

slope  WgtAvg (phase2pi , amplitude2pi , NN,  SNRO ,  NoiseFlag)  ; 


VLQx=l . /Magnitudes; 

VLQy=VLQx; 
h=l  ; 

dx=size (phase2pi, 1) /NN; 

[E, phase_CER, VLQ] =NVWCER_SHWFS_V2 (Fx, Fy, VLQx, VLQy) ; 
phase  zonal=zonal  2(Fx,Fy,h); 

[PhaseModal, Terms] =ZernikeModal (Fx*dx, Fy*dx) ; 


%%  Strehl  Ratio  Calculations  method  1; 
ims=size (phase2pi, 1) ; 

dd=ims/NN;  %Spacing  for  each  subsample  distance 


undersamp  sw=amplitude2pi (dd/ 2 : dd : ims-dd/ 2 , dd/ 2 : dd : ims- 
dd/2) . *exp (i*phase2pi (dd/2 :dd: ims-dd/2, dd/2 :dd: ims-dd/2) ) ; 
undersamp  f r=amplitude2pi ( [ 1 : dd : ims  ims ] , [ 1 : dd ; ims 
ims] ) . *exp (i*phase2pi ( [1 :dd: ims  ims] , [1 :dd: ims  ims] ) ) ; 

ref f_sw= ( sum ( sum (abs ( under samp_sw) ) ) ) A2 ; 
reff  fr= (sum (sum (abs (undersamp  fr))))A2; 
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comp  zonal=exp (~i*phase  zonal); 
comp  modal=exp (~i*PhaseModal) ; 
comp  CER  unwrap=exp (-i*phase  CER)  ; 
comp__CER_wrap=con j  (E)  ./abs  (E)  ; 

strehl  zonal= (abs (sum (sum (undersamp  sw.*comp  zonal )))) A2 /ref f  sw; 
strehl  modal= (abs (sum (sum (undersamp  sw.*comp  modal )))) A2 /ref f  sw; 
strehl  CER  unwrap= (abs (sum (sum (undersamp  fr.*comp  CER  unwrap) ))) A2 /ref f 
_fr; 

strehl  CER  wrap= (abs (sum (sum (undersamp  fr.*comp  CER  wrap) ))) A2 /ref f  fr; 

%%  Plots 

eval  (  [ 'cd  567V  num2str  (SNRO)  'V  num2str  (NN)  ]  ) 
figure; imagesc (phase  CER) ; title ([ 'CER  Unwrapped  '  'Strehl=' 
num2str (strehl  CER  unwrap)  '  SNR= '  num2str (SNRO)  'N=' 
num2str (NN) ] ) ; colorbar;  colormap  gray; 
saveas (gcf , ' CER_Unwrapped .fig', 'fig'); 

figure; imagesc (phase  zonal) ; title ([ 'Zonal  Unwrapped  '  'Strehl=' 
num2str (strehl  zonal)  '  SNR= '  num2str (SNRO)  'N=' 
num2str (NN) ]); colorbar;  colormap  gray; 
saveas (gcf, '  Zonal  Unwrapped .fig' , 'fig' ) ; 

figure; imagesc (PhaseModal) ; title ([ 'Modal  Unwrapped  '  'Strehl=' 
num2str (strehl  modal)  '  SNR= '  num2str (SNRO)  'N=' 
num2str (NN) ]); colorbar;  colormap  gray; 
saveas (gcf, ' Modal_Unwrapped .fig' , 'fig'  ) ; 

figure; imagesc (mod (phase  CER, 2*pi) ); title ([ 'CER  Wrapped  '  'Strehl=' 
num2str (strehl  CER  unwrap)  '  SNR= '  num2str (SNRO)  'N=' 
num2str (NN) ]); colorbar;  colormap  gray; 
saveas (gcf, '  CER_Wrapped .fig', 'fig'); 

figure; imagesc (mod (phase  zonal , 2 *pi )); title ([' Zonal  Wrapped  ' 

'Strehl='  num2str (strehl  zonal)  '  SNR= '  num2str (SNRO)  'N=' 
num2str (NN) ]); colorbar;  colormap  gray; 
saveas (gcf, '  Zonal  Wrapped .fig' , 'fig' ) ; 

figure; imagesc (mod (PhaseModal, 2*pi) ); title ([ 'Modal  Wrapped  '  'Strehl=' 
num2str (strehl  modal)  '  SNR= '  num2str (SNRO)  'N=' 
num2str (NN) ]); colorbar;  colormap  gray; 
saveas (gcf, ' Modal_Wrapped . fig',  'fig'); 

cd  ../../../; 


B.  SHACK-HARTMANN  WAVEFRONT  SENSOR  SIMULATION  CODE 

function  [Fx, Fy, Magnitudes, SNR]  = 

slope  WgtAvg (phase2pi , amplitude2pi , NN, SNRO , flag) 

%Create  Complex  Matrix  from  amplitude  and  phse 
Phi  Complex=amplitude2pi . *exp (li*phase2pi) ; 

%Generate  complex  differentials 

PhiDx=Phi  Complex (:, 2 : end) . *conj ( Phi  Complex (:, 1 : end-1 )) ; 

PhiDy=Phi  Complex (2  tend,  :). *conj  (Phi  Complex ( 1 : end-1 ,:)) ; 

%Magnitude  of  differentials  for  weighting 
MagX=abs (PhiDx) ; 

MagY=abs (PhiDy) ; 
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%Phase  of  differences 
thetaX=angle (PhiDx) ; 
thetaY=angle (PhiDy) ; 

Intensity=abs ( Phi  Complex). A2; 

Intensity= [ Intensity (:,  1 )  Intensity  Intensity (:, end) ] ; 

Intensity= [ Intensity ( 1 ;  Intensity;  Intensity (end, :)] ; 

%repeat  edge  values 

MagX= [MagX ( : , 1 )  MagX  MagX ( : , end) ] ; 

MagY= [MagY (1, : ) ;  MagY;  MagY (end, :)]; 
thetaX= [thetaX ( : , 1 )  thetaX  thetaX ( : , end) ] ; 
thetaY= [thetaY (1, : ) ;  thetaY;  thetaY (end, : ) ] ; 

%Calculate  differential  size 
[y,  x] =size (phase2pi) ; 
dy=round (y/NN) ; 
dx=round (x/NN) ; 

%Slope  outputs 
Fx=zeros (NN) ; 

Fy=zeros (NN) ; 

Magnitudes=zeros (NN) ; 

%Weighting  Factors 

%Since  the  column  (row)  for  x  (y)  differentials  has  overlapping  of 
sample 

%area  for  each  lenslet  and  is  ideally  located  in  the  middle  of  each 
phase 

%point  (fried  geometry)  a  weight  of  one  half  is  applied  to  the  over 
sampled 

%values  on  the  edges 
Wx=ones (dy, dx+1 ) ; 

Wx ( : , 1 ) =0 . 5 ; 

W x ( : , end) =0.5; 

Wy=ones (dy+1 , dx) ; 

Wy (1, : ) =0 . 5; 

Wy (end, : ) =0 . 5 ; 

Wmag=ones (dy+1 , dx+1 ) ; 

Wmag ( : , 1) =0 . 5 ; 

Wmag ( : , end) =0.5; 

Wmag ( 1 , : ) =Wmag ( 1 , : ) . * . 5 ; 

Wmag (end, : ) =Wmag (end, : ) . * . 5 ; 

thetaFx=zeros (dy,  dx+1 ) ; 
thetaFy=zeros (dy+l,dx)  ; 

magFx=zeros (dy, dx+1 ) ; 
magFy=zeros (dy+1 , dx) ; 
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for  u— 1 : NN 

for  v=l : NN 

thetaFx=thetaX ( (v-1) *dy+l :v*dy, (u-1) *dx+l :u*dx+l) ; 
thetaFy=thetaY ( (v-1) *dy+l :v*dy+l, (u-1) *dx+l :u*dx) ; 
magFx=MagX ( (v-1) *dy+l :v*dy, (u-1) *dx+l :u*dx+l) . *Wx; 
magFy=MagY ( (v-1) *dy+l :v*dy+l, (u-1) *dx+l :u*dx) . *Wy; 

Fx (v, u) =sum ( sum ( thetaFx . *magFx) ) / sum ( sum (magFx) ) . *dx; 
Fy (v, u) =sum (sum (thetaFy . *magFy) ) / sum (sum (magFy) ) . *dy; 
Magnitudes (v, u) =sum ( sum ( Intensity ( (v-1 ) *dy+l : v*dy+l , ( 
1) *dx+l :u*dx+l) . *Wmag) ) ; 
end 

end 

SNR=SNRO*Magnitudes/mean (Magnitudes ( : ) ) ; 

delX=randn (NN, NN) . *pi . /SNR; 
delY=randn (NN, NN) . *pi . /SNR; 

if  flag==l 

Fx=Fx+delX; 

Fy=Fy+delY; 

end 

end 


C.  ZONAL  RECONSTRUCTION 

function  [phases]  =  zonal  2(Fx,Fy,h) 

[a,  b] =size (Fx) ; 

[c,  d] =size (Fy) ; 

if  c~=a  | |  d~=a 

error ('Size  of  Fx  differs  from  Fy' ) ; 

end 

if  a~=b 

error ('Fx  is  not  square'); 

end 

if  c~=d 

error ('Fy  is  not  square'); 

end 

N=a; 

S= [reshape (Fx'  ,  1 ,  NA2 )  reshape (Fy'  ,  1 ,  NA2 ) ]  '  ; 

A=f ormA (N) ; 

D=f ormD (N) ; 

phase s=pinv (A) *D*S ; 

%  [U,  A2,  V]  =svd  (A,  0)  ; 
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%tiltx 

%tilty 


A2=pinv (A2 ) ; 
phases=V*A2  *U' *D*S; 


phase s=re shape (phases , N, N) '  . /h; 
end 

%  Derived  from  (Dai  2008) 
function  A=formA(N) 

A=zeros (2*N* (N-l) ,NA2) ; 
for  i=l : N 

for  j=l : (N-l) 

A  (  ( i — 1 ) * (N-l) +j,  (i-1) *N+j ) =-l ; 

A  (  ( i-1 ) * (N-l ) + j ,  (i-1) *N+j+l)=l; 

A ( (N+i-1 ) * (N-l) +j , i+ (j-1) *N)=-1; 

A (  (N+i-1) * (N-l) +j , i  +  j  *N) =1; 

end 

end 

end 

function  D=formD(N) 

D=zeros (2*N* (N-l) , 2*NA2) ; 
for  i=l : N 

for  j=l : (N-l) 

D(  (i-1) * (N-l) +j ,  (i-1) *N+ j ) =0 . 5 ; 

D  (  (i-1) * (N-l) +j ,  (i-1) *N+j+l) =0.5; 

D  (  (N+i-1) * (N-l) +j ,N* (N+j-1) +i)  =0 . 5; 
D ( (N+i-1) * (N-l ) + j , N* (N+j ) +i) =0.5; 

end 

end 

end 


D.  MODAL  RECONSTRUCTION 

function  [ PhasesModal , terms] =ZernikeModal (Fx, Fy) 


9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9-9'9'9'9-9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'S- 

oooooooooooooooooooooooooooooooooooooooooooooooo 

9-  9-  9-  9- 
o  o  o  o 

%%  Adapted  from  Matthew  Allen  Modal  Code  (c) 

?nn7  2'2'2'2'2'2'2'&-2'2'2'2'2'2'2'2'2'2'2'2'2'2'2'2'2'2' 

tUu  /  oooooooooooooooooooooooooo 


9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9-9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9' 

oooooooooooooooooooooooooooooooooooooooooooooooo 

9-  9-  9-  9- 
o  o  o  o 


NN=size (Fx, 1 ) ; 

[x  ,  y] =meshgrid ( -1+ (2 / (2*NN) )  : 2 / (NN)  : 1 )  ; 
x=reshape (x,  [  ] , 1 ) ; 
y=reshape (y,  [],!); 


S= [reshape (Fx,  [  ] , 1) ;  reshape (Fy,  [],!)]; 
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X 

2 

=  X  . 

CM 

< 

X 

*3 

=  X  . 

A3  ; 

X 

'4 

=  X  . 

A4 ; 

X 

"5 

=  X  . 

LO 

< 

y_ 

'2 

=  y- 

CM 

< 

y_ 

"3 

=  y- 

A3; 

y_ 

'4 

=  y- 

A4 ; 

y_ 

"5 

=  y* 

LO 

< 

dZdx=zeros (size (x,  1)  ,21) ; 
dZdy=zeros (size (y, 1)  ,21)  ; 

%Partial  Derivatives  Evaluation  dX 

dZdx ( : , 1) =1 ; 

dZdx ( : , 2 ) =0 ; 

dZdx ( : , 3 ) =4 . *x; 

dZdx ( : , 4 ) =2 . *x; 

dZdx ( : , 5) =2 . *y; 

dZdx ( : , 6) =-2+9 . *x  2+3. *y  2; 

dZdx (:,7)=6.*x.*y; 

dZdx ( : , 8) =-12 . *x+24 . *x. * (x  2+y  2); 
dZdx ( :  ,  9) =3 . *x_2-3 . *y_2 ; 
dZdx ( : , 10 ) =6 . *x . *y; 

dZdx ( : , 1 1 ) =-6 . *x+8 . *x . * (x  2+y  2)+8.*x  3-8.*x.*y  2; 
dZdx ( : , 12 ) =-6 . *y+8 . *y . * (x  2+y  2)+16.*x  2.*y; 

dZdx ( : , 13) =3-36 . *x_2-12 . *y_2+10 . * (x_2+y_2 ) . A2+40 . *x_2 . * (x_2+y_2 ) 
dZdx ( : , 14) =-24 . *x. *y+40 . *x. *y. * (x  2+y  2); 
dZdx(:,15)=24.*x-120.*x.* (x_2+y_27+120 . *x . * (x_2+y_2) .  A2; 
dZdx ( : , 16) =4 . *x_3-12 . *x . *y_2 ; 
dZdx ( : , 17) =12 . *x_2 . *y-4 . *y_3 ; 

dZdx ( : , 18) =-12 . *x_2+12 . *y_2+15 . *x_2 . * (x_2+y_2 ) +10 . *x_4- 
15.*y_2.*(x_2+y_2)-30.*x_2.*y_2; 

dZdx ( :  ,  19) =-24 . *x. *y+30 . *x. *y. * (x  2+y  2)+30.*x  3 . *y-10 . *x. *y  3; 
dZdx(:,20)=12.*x-40.*x.* (x_2+y_2 ) - 

40.*  x_3  +  4  0 . *  x . *  y_2  +  3  0 . *  x . * (x_2+y_2)  . A2  +  60 . *x_3 . * (x_2+y_2) - 

60.*x.*y  2.*(x  2+y  2); 

dZdx ( : ,21) =12 . *y-40.*y.* (x_2+y_2) - 

80.*x_2.*y+30.*y.* (x_2+y_2) . A2+120*x_2 . *y . * (x_2+y_2) ; 

%  dZdx ( sqrt (x_2+y_2 ) >1 ,  : ) =0  ; 

%Partial  Derivatives  Evaluation  dY 

dZdy ( : , 1 ) =0 ; 

dZdy ( : , 2) =1; 

dZdy ( : , 3) =4 . *y; 

dZdy ( : , 4 ) =-2 . *y; 

dZdy ( : , 5) =2 . *x; 

dZdy ( : , 6) =6 . *x . *y; 

dZdy ( :  ,  7 ) =-2  +  3 . *x_2  +  9 . *y_2 ; 

dZdy ( : , 8) =-12 . *y+24 . *y . * (x_2+y_2) ; 

dZdy ( : , 9) =-6 . *x . *y; 

dZdy ( : , 10) =3 . *x_2-3 . *y_2; 

dZdy ( : , 1 1 ) =6 . *y+8 . *x  2 . *y-8 . *y . * (x  2+y  2)-8.*y  3; 
dZdy ( : , 12 ) =-6 . *x+8 . *x . * (x  2+y  2)+16.*x.*y  2; 
dZdy ( : , 13) =-24 . *x. *y+40 . *x. *y. * (x  2+y  2); 

dZdy ( : , 14 ) =3-12 . *x_2-36 . *y_2+10 . *7x_2+y_2 ) . A2+4 0 . *y_2 . * (x_2+y_2) 


dZdy ( : , 15) =2  4 ,*y-120.*y.* (x_2+y_2) +120 . *y. * (x_2+y_2) ,A2; 

dZdy(:,16)=-12.*x_2.*y+4.*y_3; 

dZdy ( : , 17 ) =4 . *x_3-12 . *x . *y_2 ; 

dZdy ( : , 18) =24 . *x. *y+10 . *x_3 . *y-30 . *x. *y. * (x  2+y  2)-30.*x.*y  3; 
dZdy ( : , 19) =-12 ,*x_2+12 . *y~2+15 . *x_2 . * (x_2+y~2 ) +30 . *x_2 . *y_2r 
15.*  y_2 . * ( x_2 +y_2 ) - 1 0 . *  y_4 ; 
dZdy ( : , 20) =-12 . *y- 

40. *x  2 . *y+4 0 . *y . * (x  2+y_2)+40.*y  3+60. *x  2.*y.*(x  2+y  2)- 
30.*y.*(x  2+y  2 ) . A2-60 . *y_3 . * (x  2+y  2); 
dZdy ( : , 2l7=127*x-40 . *x. * (x_2+y_2) - 

8  0 . *  x . *  y_2  +  3  0 . *  x . * (x_2+y_2)  . A2  +  120 . *x . *y_2 . * (x_2+y_2) ; 

%  dZdy ( sqrt (x_2+y_2 ) >1 ,  : ) =0 ; 

dZ= [dZdx; dZdy] ; 

terms=pinv (dZ) *S; 
a=terms ; 

[x,  y]=meshgrid( (-1+2/NN/2) :2/NN: (1-2/NN/2) ) ; 

%  [x  y] =meshgrid ( -1 : 2 /NN : 1 ) ; 


X 

2 

=  X  . 

CN 

< 

X 

3 

=  X  . 

A3  ; 

X 

A 

=  X  . 

A4 ; 

X 

*5 

=  X  . 

LO 

< 

y_ 

'2 

=  y- 

CM 

< 

y_ 

"3 

=  y- 

A3  ; 

y_ 

4 

=  y- 

A4 ; 

y_ 

*5 

=  y- 

LO 

< 

PhasesModal=a ( 1 ) *x+ . .  . 
a  (2 ) *y+ . . . 

a (3) * (-1+2 .* (x_2+y_2) )  +  .  . . 
a ( 4 ) * ( x_2 - y_2 ) + .  .  . 
a (5) * (2 . *x. *y) + . . . 
a (6) * (-2 . *x+3 . *x. * (x  2+y  2))+... 
a ( 7  )  * ( -2 . *y+3 . *y . * (x  2+y  2))  +  ... 
a (8) * (1-6.* (x_2+y_2) +6. *Tx_2+y_2)  . A2)  +  .  . . 
a (9) * (x_3-3.*x.*y. A2) +. . . 
a  (10) * (3. *x_2 . *y-y_3)  +  . . . 

a  (11) * (-3. *x_2  +  3.*y_2  +  4 ,*x_2 .* (x_2+y_2) -4 ,*y_2 .* (x_2+y_2) ) +. . . 

a  ( 12 ) *  ( -6 . *x . *y+8 . *x . *y . * (x  2+y  2))  +  ... 

a (13) * (3 . *x-12 . *x. * (x_2+y_2T+10T*x.* (x_2+y_2) . A2) +. . . 

a (14) * (3.*y-12 ,*y.* (x_2+y_2) +10 . *y. * (x_2+y_2) . A2) +. . . 

a (15) * (-1+12 . * (x_2+y_2 ) -30.* (x_2+y_2) . A2+20.* (x_2+y_2) . A3) +. . . 

a (16) * (x_4-6.*x_2 ,*y_2+y_4) +. . . 

a (17) * (4 ,*x_3.*y-4 ,*x.*y_3) +. . . 

a  (18)  *  (~4.*x_3  +  12.*x.*y_2  +  5.*x_3.*  (x__2+y_2)  -15  .  *x  .  *y_2  .  *  (x_2+y 
a (19) * (-12 . *x  2.*y+4.*y  3+15. *x  2.*y.*(x  2+y  2 ) -5 . *y_3 . * (x  2+y 
a (20 )  * ( 6 . *x_2-6 . *y_2- 

20 . *x_2 . * (x_2+y_2 ) +20 . *y_2 . * (x_2+y_2 ) +15 . *x_2 . * (x_2+y_2) ,A2- 
15 . *y_2 . * (x_2+y_2 )  . A2 )  +  .  .  . 
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CM  CN 


a(21) * (12.*x.*y-40.*x.*y.* (x_2+y_2 ) +30 . *x . *y . * (x_2+y_2) .  A2)  ; 

%  PhasesModal ( sqrt (x_2+y_2 ) >1 ) =0 ; 
end 


E.  NVWCER 

function  [E, phi , VLQ] =NVWCER  SHWFS  V2 (dpx, dpy, sigsqx, sigsqy, mask, amp) 

S-S-S'S-S'S-S'S-S'S-S'S-S'S-S'S-S'S-S'S-S-S-S'S-S-S-S'S-S'S-S'S-S'S-S'S-S'S-S'S-S'S-S'S-S'S-S'S-S'S-S'S-S'S-S'S-S'S-S'S-S'S-S'S-S'S-S'S- 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

0,0,  0,0, 
o  o  o  o 

%%%%%%  Adapted  from  original  code  provided  by  Dr.  David  Fried 

9-9-9'9-9'9'9'9'9'9- 

oooooooooo 

%%%%%%  All  credit  for  original  code  is  given  to  him  and  provided 

TT-i  -HViS-S-S-S-S-S- 
W1  LII  o  o  o  o  o  o 

%%%%%  his  knowledge 

9'9'9'9'9'9'9'9'9'9-9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9-9'9'9'9'9'9'9'9-9'9-9'9'9'9'9'9-9'9'9'9-9'9- 

oooooooooooooooooooooooooooooooooooooooooooooooooooooo 


%%  Error  Checking  for  Size  of  input  matricies 
[a,  b] =size (dpx) ; 
if  a~=b 

error ('dpx  is  not  a  square  array.'); 

end 

%  N=round (log2 (a) ) ; 

%  if  a~=2 AN 

%  error ('Size  of  dpx  is  not  2 AN-by-2 AN . ' ) ; 

%  end 

[c, d] =size (dpy) ; 
if  c~=a  |  d~=a 

error ('Size  of  dpy  is  different  from  that  of  dpx . ' ) ; 

end 

[c,d]=size (sigsqx)  ; 
if  c~=a  |  d~=a 

error ('Size  of  sigsqx  is  different  from  that  of  dpx.'); 

end 

[c,d]=size (sigsqy)  ; 
if  c~=a  |  d~=a 

error ('Size  of  sigsqy  is  different  from  that  of  dpy.'); 

end 

if  nargin<4 

error ( 'Not  enough  inputs .  You  atleast  need  X&Y  phasers  and  X&Y 
noise' ) ; 
end 

if  nargin==4 

mask=ones (a+1)  ; 

end 
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%%%  Assign  Hartman  data  to  first  lattace  and  carry  out  reconstruction 

o  o  o  o  o 


9'9'9'9'5'9'5'9'9'9'5'9'2'9'5'9'2'9'5'9'2'9'9'9'9'9'5'9'9'9'5'9'9'9'9'9'2'9'2'9'9'9'9'9'2'9'2'9'2'9'5'S-2'9'5'9'2'9'9'9'9'9'5'S-9'9'9'9'9'9'2' 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

0.00  O. 

o  o  o  o 

9'9'9'9'9'9-9'9'9'9'9'S-9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'S-9'9'9'9-9'9'9'9'9'9'9'9-9'9'9'9'9'9'9'9'9'9'9'9-9'9'9'&-9'9'9'9-9'9'9'S-9'9'9'9'9'9'9' 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


NN=ceil (log2 (a) ) ;  %Creates  propper  size  2AN  lattice 

L=2ANN;  %Number  of  unit  squares  to  be  used 

dpu=zeros (L+l , L) ; 

dpv=zeros (L, L+l ) ; 

sigsqu=inf *ones (L+l ,  L)  ; 

sigsqv=inf *ones (L, L+l )  ; 

%assigning  values  from  Shack  Hartmann  Sensor  to  first  lattice  for 
%exponential  reconstructor, 
for  x=l : a 

for  y=l:a 

if  x+y==2 *round ( (x+y) /2 ) 
u— (x-y+L) /2+1 ; 
v= (x+y) / 2 ; 

dpv (v, u) =dpx (y, x) +dpy (y, x) ; 
sigsqv(v,u)=sigsqx(y,x)+sigsqy(y,x); 

else 

u— (x-y+L+1 ) /2 ; 
v— (x+y+1 ) / 2  ; 

dpu (v, u) =dpx (y, x) -dpy (y,  x) ; 
sigsqu(v,u)=sigsqx(y,x)+sigsqy(y,x); 

end 

end 

end 

%Hudgin  lattice  generated  for  phase  1.  Passed  into  NWRCER  with 
differential 

%phasors  calculated  in  the  equation  (exp(i*dp#)) 

[Ep, sigsqp] =NVWCER4SID (exp (i*dpu) , exp (i*dpv) , sigsqu, sigsqv) ; 

%Create  raw  lattices  so  that  only  values  which  correpsond  to  the 
orignal 

%lattice  are  taken  from  the  hudgin  geometry  generation  and  calculation 
El=NaN*ones (a+1 , a+1 ) ; 
sigsql=NaN*ones (a+1 ,  a+1 )  ; 

%Go  through  point-by-point  on  large  (hudgin/phase  1  lattice)  and  assign 
%only  values  which  correspond  to  values  overlapping  with  the  original 
%lattice.  If  they  are  outside  the  appearture  they  maintain  the  value  of 
0 

%phase  and  Inf  variance  (noise) 
for  u— 1 : L+l 
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for  v=l : L+l 

x=u+v-l-L/ 2 ; 
y=-u+v+l+L/ 2 ; 

if  ( (x>=l  &  x<=a+l)  &  (y>=l  &  y<=a+l))  &  isf inite ( sigsqp (u, v) ) 

El  (y,  x)  =Ep  (v,  u)  ; 
sigsql (y, x) =sigsqp (v, u) ; 

end 

end 

end 


9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'&-9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'&-9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'S-9'9'9'9'9'9'9'9'9'9'9'9-9'9' 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

0,0,  0,0, 
o  o  o  o 

9'9-9'9-9'9'9'9'9'9-9'9-9'9'9'9'9'9-9'9'9'9'9'9-9'9-9'9-9'9'9'9'9'9-9'9-9'9'9'9'9'9-9'S-9'9-9'9-9'9-9'9'9'9'9'9-9'9-9'9'9'9-9'9'9'9'9'9-9'9- 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


0,0,  0,0, 
o  o  o  o 

%%%  Assign  Hartman  data  to  second  lattace  and  carry  out  reconstruction 

0,0,  0,0, 
o  o  o  o 

9'9'9'9'5'9'5'9'9'9'5'9'5'9'9'9'9'9'9'9'5'9'2'9'2'9'9'9'5'9'5'9'9'9'9'9'5'9'5'9'5'9'2'9'5'S-9'9'9'9'9'9'2'S-2'9'9'9'9'9'2'9'2'9'9'9'9'9'9'9' 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

0,0,  0,0, 
o  o  o  o 


9'9'9'9'9'9'9'S-9'9-9'9'9'9-9'9'9'9'9'9'9'9-9'S-9'9'9'9'9'9'9'9'9'9'9'9'9'9-9'S-9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9-9'9'9'9'9'9-9'9'9'9' 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

0,0,  0,0, 
o  o  o  o 


dpu=zeros (L+l ,  L)  ; 
dpv=zeros (L, L+l ) ; 
sigsqu=inf *ones (L+l ,  L)  ; 
sigsqv=inf *ones (L, L+l )  ; 

for  x=l : a 

for  y=l:a 

if  x+y==2 *round ( (x+y) /2 ) 
u— (x-y+L) / 2 ; 
v= (x+y) / 2 ; 

dpu (v, u) =dpx (y, x) -dpy (y, x) ; 

sigsqu (v, u) =sigsqx (y, x) +sigsqy (y, x) ; 

else 

u— (x-y+L+1 ) /2 ; 
v= (x+y-1 ) / 2 ; 

dpv  (v,  u)  =dpx  (y,  x)  +dpy  (y,  x)  ; 
sigsqv(v,u)=sigsqx(y,x)+sigsqy(y,x); 

end 

end 

end 

[Epp, sigsqpp] =NVWCER4SID (exp (i*dpu) , exp (i*dpv) , sigsqu, sigsqv) ; 

E2=NaN*ones (a+1 , a+1 ) ; 
sigsq2=NaN*ones (a+1 ,  a+1 )  ; 

for  u— 1 : L+l 

for  v=l : L+l 

x=u+v-L/ 2 ; 
y=-u+v+l+L/ 2 ; 

if  ( (x>=l  &  x<=a+l)  &  (y>=l  &  y<=a+l))  &  isf inite ( sigsqp (u, v) ) 

E2  (y,  x)  =Epp  (v,  u)  ; 
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sigsq2 (y,  x) =sigsqpp (v, u) ; 

end 

end 

end 


9'9'9'9'5'9'5'9'9'9'5'9'9'9'9'9'2'9'2'S-9'S-2'9'9'9'9'9'5'9'9'9'2'9'5'9'5'9'5'9'5'9'9'9'9'9'9'9'2'9'9'9'5'9'9'9'5'9'9'9'5'9'9'9'2'9'2'S-9'9'2' 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

0,0.  0,0, 
o  o  o  o 

S-9'9'9'2'9'9'9'9'9'9'9'5'9'2'9'9'9'9'9'5'9'5'9'9'9'5'9'2'9'5'9'9'9'9'9'5'9'2'9'2'9'9'9'5'9'5'9'9'9'5'9'5'9'2'9'9'9'9'9'9'9'5'9'9'9'9'9'5'9'9' 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

0,0,  0,0, 
o  o  o  o 

%%%  MERGE  MATRICES  BACK  TOGETHER  INTO  ORGINAL 

9'9'9'9'9'S-9'9'9-9-9-9-S-9'9'9'9'9'9'9'9-9'9'9'9'9'9' 

ooooooooooooooooooooooooooo 

9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9-9'9'9'9'9'9-9'S-9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9-9'9'9'9-9'9'9'9'9'9-9' 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

0,0,  0,0, 
o  o  o  o 

9-9-9'9'9'9'9'9-9'9'9'9-9'9'9'S-9'9'9'9-9'9'9'9'9'9-9'9'9'9'9'9-9'9-9'9-9'9'9'9-9'9-9'9-9'9'9'9-9'9'9'9'9'9'9'9-9'9-9'9-9'9'9'S-9'9-9'9'9'9-9' 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

0,0,  0,0, 
o  o  o  o 


E2t=NaN*ones (a+1) ; 
sigsq2t=NaN*ones (a+1)  ; 
for  x=l : a+1 

for  y=l : a+1 

if  x+y~=2 *round ( (x+y) /2  ) 

S=0 ;  T=0 ; 
if  x— 1>— 1 

S=E1 (y, x-1) /sigsql (y, x-1) ; 

T=1 / sigsql (y,x-l) ; 

end 

if  x+l<=a+l 

S=S+E1 (y, x+1) /sigsql (y, x+1) ; 

T=T+1/ sigsql (y, x+1) ; 

end 

if  y-l>=l 

S=S+E1 (y-l,x) /sigsql (y-l,x) ; 

T=T+l/sigsql (y-l,x) ; 

end 

if  y+l<=a+l 

S=S+E1 (y+l,x) /sigsql (y+l,x) ; 

T=T+l/sigsql (y+l,x) ; 

end 
T— 4 / T ; 

S=S/ (4/T) ; 

E2t (y, x) =S ; 
sigsq2t (y, x) =T; 

end 

end 

end 

[x, y] =meshgrid ( 1 : a+1 ) ; 

ii=f ind (x+y~=2*round ( (x+y) / 2) ) ; 

S=sum (E2t(ii) .*conj (E2 (ii) ) . *mask (ii) ./ (sigsq2t(ii)+sigsq2 (ii) ) ) ; 
phi21=angle (S) ; 

E=E  1  ; 
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E  ( ii ) =E2 ( ii ) *exp ( i*phi2 1 ) ; 
VLQ=sigsql ; 

VLQ (ii) =sigsq2 (ii) ; 
jj=find(~isfinite (VLQ) ) ; 

E ( j j ) =NaN; 


E=E . /abs (E) ; 
if  nargin==6 
E=amp . *E ; 

end 

phi=BranchPointPhase (E)  ; 
end 


9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9-9'9'9'9'9'9'9'9'9'9'9'9-9'9-9'9'9'9'9'9'9'9-9'9'9'9'9' 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

S'9-9'9'9- 
o  o  o  o  o 

S-  9- 
o  o 

%  Function,  called  BranchPointPhase,  used  to  generate  a  phase  function, 
PHI, 

%  that  corresponds  (modulo  2*pi)  to  the  phase  of  the  input  complex 
phasor , U, 

%  and  that  is - as  far  as  possible - smooth,  i.e.,  which  minimizes 

difference 

%  between  values  of  phi  between  adjacent  points  in  the  array  of  values 
of 


%  phi.  Ideally,  the  magnitude  of  the  difference  of  phase  values  between 
%  adjacent  points  should  be  less  than  pi  everywhere  except  at  branch- 
cuts  . 


%  This  function  is  formed  as  the  combination  of  the  two  functions 
SmoothPhase 

%  and  Reconstructor  that  are  listed  in  TN-100.  The  contribution  that 
comes 

%  from  (or  rather  that  was)  the  function  called  Reconstructor  appears 
here  as 

%  a  sub  function,  clearly  identified  as  being  Reconstructor,  in  the 
final 

%  part  of  this  listing.  All  the  rest  of  the  listing  comes,  essentially 
%  intact,  from  SmoothPhase. 


%  The  input  function,  U,  defines  the  extent  of  the  aperture  by  its 
value 

%  being  equal  to  NaN  for  those  array  points  corresponding  to  positions 
%  outside  the  aperture.  The  array  on  which  U  is  defined  is  to  be 
square,  with 

%  U  corresponding  to  a  clear  circular  aperture  just  smaller  than  the 
size  of 

%  the  array.  For  points  on  the  array  that  are  within  the  aperture  the 
values 

%  of  U  are  complex  numbers  (representing  phasors  of  the  optical  field) 
while 

%  for  points  on  the  array  that  are  outside  the  aperture  the  values  of  U 


are 
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NaN. 


O, 

o 

%  The  calculation  of  PHI  is  accomplished  by  first  evaluating  the 
difference 

%  of  the  phase  (modulo  2*pi)  between  adjacent  values  of  U,  using  those 
phase 

%  differences,  dPx  and  dPy,  to  calculate  the  curl  for  each  elemental 
square 

%  of  the  array,  and  using  the  values  of  the  curl  to  determine  the 
location 

%  and  sign  of  the  branch  points.  Using  this  branch  point  information 
the 

%  complex  phaseor,  Uh,  associated  with  the  branch  points  is  calculated 
and 

%  from  this  the  hidden  phase,  phi,  is  then  calculated.  The  ratio  U  over 
Uh, 

%  denoted  by  u,  is  calculated.  The  phase  differences  dpx  and  dpy 
associated 

%  with  u  are  evaluated.  Then  the  associated  phase  is  generated  by  first 
%  adding/subtracting  dpy  values  along  the  centeral  line  of  the  array  to 
give 

%  phase  values  for  the  all  the  points  above/below  the  center  point  of 
the 

%  central  line,  and  then  starting  from  each  point  on  that  central  line 
by 

%  adding/subtracting  dpx  values  the  phase  value  for  each  point  to  the 
%  right/left  of  the  central  point  line  the  phase  value  for  each  point 
in  the 

%  array  is  calculated.  To  these  phase  values,  which  may  be  associated 
with 

%  the  non  hidden  phase  values,  the  hidden  phase  values,  phi,  are  added 
to 

%  yield  the  final  output  phase  values,  PHI. 

o, 

o 

%  For  those  cases  in  which  there  the  number  of  positive  sign  branch 
points  is 

%  not  equal  to  the  number  of  negative  sign  branch  points  (presumably 
because 

%  the  "missing"  branch  point  is  outside  the  aperture  and  so  is  not 
observed) 

%  an  additional  branch  point  of  the  appropriate  sign  is  provided.  This 
branch 

%  point  is  placed  at  a  position  just  outside  the  aperture  and  at  a 
location 

%  chosen  on  the  basis  of  a  "potential  field, "  V,  formed  by  the  positive 
and 

%  negative  sign  branch  points.  The  position  is  chosen  to  correspond  to 
the 

%  maximum  of  the  potential  field - the  maximum  over  positions  just 

outside 

%  the  aperture.  This  added  branch  point  is  introduced  to  "guide"  the 
branch 

%  cut  that  goes  out  of  the  aperture. 

9'9'9'9'9'9-9'9-9'9'9'9'9'S-9'9-9'9'9'9'9'9'9'9-9'9'9'9'9'9'9'9'9'9-9'9'9'9'9'9'9'9'9'9'9'9'9'9-9'9'9'9'9'9'9'9-9'9'9'9'9'9'9'9'9'9'9'9'9'9'9' 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

g,  9-  9-  9-  9- 
o  o  o  o  o 
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NaN 

O, 

o 

the 


%  PHI 
to  the 


INPUTS 

(N,M)  =  Complex  phasor  array.  (The  values  of  U  are  equal  to 
for  those  array  points  whose  position  is  outside 
aperture . ) 

OUTPUTS 

(N,M)  =  Maximally  smooth  phase  corresponding  (modulo  2*pi) 
phase  of  U.  (Values  of  PHI  are  equal  to  NaN  for 


points 


outside  the  aperture.) 

BPcount  (1,1)  =  Number  of  branch  points  in  U  (located  within  the 
aperture) .  respectively. 

(N,M)  =  Hidden  phase  corresponding  to  U.  (Values  of  phi  are 


%  phi 
equal 

o, 

o 

%  BPes 
column 


(i.e. 

O, 

o 

third 

O, 

o 

the 

o, 

o 

a  is 

o, 

o 

those 


to  NaN  for  points  outside  the  aperture.) 

(a, 3)  =  The  second  column  of  values  gives  the  x-axis  (i.e., 

number)  position  of  a  brach  point,  while  the 
corresponding  first  column  values  give  the  y-axis 

the  row  number)  position  of  the  branch  point.  The 

column  value  is  equal  to  +1  or  to  -1,  indicating  if 

branch  point  is  positive  or  negative.  (The  value  of 

equal  to  the  number  of  branch  points,  including 

added  at  positions  outside  the  aperture.) 

[PHI, BPcount, phi, BPes] =BranchPoint Phase (U) 


function  [ PHI , BPcount, phi , BPes ] =BranchPoint Phase (U) 

[N,  M]  =size  (U)  ; 

if  N~=M;  error ( 'Array  should  be  square.');  end 
%Equations  35a/b 

%Atan2  is  for  arg ( )  which  calculates  principle  value  form 

dPx=U ( : , 2 : end) . *conj (U ( : , 1 : end-1 ) ) ; 

dPx=atan2 ( imag (dPx) , real (dPx) ) ; 

dPy=U (2 : end,  : )  . *conj  (U ( 1 : end-1,  : ) )  ; 

dPy=atan2 (imag (dPy) , real (dPy) ) ; 

BPes= [ ] ; 

curl=dPx ( 1 : end-1 ,  : ) +dPy ( : , 2 : end) -dPx (2 : end,  : ) -dPy ( : , 1 : end-1 ) ;  %Eqn  36 
for  n=l : N-l 

for  m=l :M-1 


80 


%  if  abs (curl (n, m) ) >le-3  %Tollerace  for  Branch  point  from 

Curl  of  unit 

if  abs (curl (n, m) )>. l*pi 

BPes=[BPes;  n  m  sign (curl (n, m) )] ; 

end 

end 

end 

BPcount=size (BPes, 1)  ; 
bpcount=BPcount ; 

Uh=ones (N, M) ; 
if  bpcount>0 

BPexcess=sum (BPes  ( : , 3 ) )  ; 

%Concerned  with  ensuring  the  number  of  positive  and  negative  branch 
%points  are  the  same.  And  if  not  carry  out  this  process 
while  BPexcess~=0  %Calculating  Step  2  1/2  (eqn  44) 

R= (N+3 ) /2 ; 

theta= (0:359) *pi/180; 
x=R*cos (theta) +M/2; 
y=R*sin (theta) +N/2; 

V=zeros (1, 360) ; 
for  k=l:BPcount 

V=V+BPes (k, 3) .  / sqrt ( (x-BPes (k, 2) ) ,A2+(y-BPes(k,l) ) ,A2) ; 

end 

[mx, ii ] =max (BPexcess*V) ; 
bpcount=bpcount+l ; 

BPes=[BPes;  y(ii)  x(ii)  -sign (BPexcess ) ] ; 

BPexcess=sum (BPes ( : , 3 ) ) ; 

end 

[x, y] =meshgrid (1 :M, 1 :N) ; 
for  bpc=l : bpcount 

X=BPes (bpc, 2) +0 . 5; 

Y=BPes (bpc, 1 ) +0.5; 
pn=BPes (bpc, 3) ; 
if  pn>0 

Uh=Uh . * [ (x-X) +i* (y-Y) ] ; 

else 

Uh=Uh . / [ (x-X) +i* (y-Y) ] ; 

end 

end 

end 

phi=angle (Uh) ; 

ii=f ind ( ~isf inite (U) ) ; 

phi ( ii ) =NaN; 


u=U. /Uh; 

dpx=u ( : , 2 : end) . *con j (u ( : , 1 : end-1 ) ) ; 
dpx=atan2 ( imag (dpx) , real (dpx) ) ; 
dpy=u (2 : end, : ) . *conj (u ( 1 : end-1, : ) ) ; 
dpy=atan2 (imag (dpy) , real (dpy) ) ; 
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Phi=Reconstructor (dpx, dpy) ; 

PHI=Phi+phi ; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
Code  provided  by  Dr.  David  Fried  as  apart  of  NVWCER  algorithm  fo 
%  Function,  called  Reconstructor,  used  to  accomplish  reconstruction 
based  on 

%  simple  addition  of  phase  differences.  Starting  from  the  center  of  the 
array 

%  first  the  phase  differences  along  the  y-axis  direction  are  added  for 
the 

%  central  line.  Then  starting  from  each  point  in  that  line  the  phase 
%  differences  are  added  along  the  x-axis  direction. 

%  The  process  starts  with  the  phase  of  the  central  element  treated  as 
being 

%  equal  to  zero,  but  finally  the  phase  of  all  elements  are  adjusted  so 
that 

%  mean  phase  is  equal  to  zero. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%% 

% 

%  INPUTS 

%  pdx  (N,N-1) 

% 

% 

%  pdy  (N-1,N) 

% 

radians) . 

% 

%  OUTPUTS 

%  phi  (N,N) 

% 

%  phi=Reconstructor (pdx,  pdy) 

% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%% 

function  phi=Reconstructor (pdx,  pdy) 

if  nargin==2 
yn=-l ; 

end 

[a, b] =size (pdx) ; 
if  a~=b+l 

error ('Size  of  pdx  is  not  N-by- (N-l) . ' ) 

end 

[c, d] =size (pdy) ; 
if  (d~=a)  I  (c~=b) 

error ('Size  of  pdy  does  not  properly  correspond  to  that  of  pdx.') 

end 

82 


Phase  difference  for  the  x-direction,  i.e.,  for  the 
direction  in  which  the  column  numbers  change  (in 
radians) . 

Phase  difference  for  the  y-direction,  i.e.,  for  the 
direction  in  which  the  row  numbers  change  (in 


Reconstructed  phase  (in  radians) . 


N=size (pdx, 1 ) ; 

hN=round (N/2 ) ; 

pd=-f lipud (pdy ( 1 : hN-1 , hN) ) ; 

phi=f lipud (cumsum (pd) ) ; 

pd=pdy (hN :N-1, hN) ; 

phi=[phi;  0;  cumsum(pd) ] ; 

pd=-f liplr (pdx ( : , 1 : hN-1 )  )  ; 
phi=fliplr (cumsum ( [phi  pd],2)); 
pd=pdx ( : , hN : N- 1 )  ; 

phi= [phi (:, 1 : hN-1 )  cumsum ( [phi (:, hN)  pd],2)] 

W=isfinite (phi) ; 
phi=phi-mean (phi (W (:))); 
phi ( ~W) =0 ; 
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